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ABSTRACT 

Motivated by the New Horizons mission, we consider how Pluto's small satel- 
lites - currently P5, Nix, P4, and Hydra - grow in debris from the giant impact 
that forms the Pluto-Charon binary or in solid material captured from the pro- 
toplanetary debris disk. If the satellites have masses close to their minimum 
masses, our analysis suggests that capture of material into a circumplanetary or 
circumbinary debris disk is a viable mechanism for satellite formation. If the 
satellites are more massive, they probably form in debris from the giant impact. 
After the impact, Pluto and Charon accrete some of the debris and eject the rest 
from the binary orbit. During the ejection, high velocity collisions among de- 
bris particles produce a collisional cascade, leading to the ejection of some debris 
from the system and enabling the remaining debris particles to find stable orbits 
around the binary. Our numerical simulations of viscous diffusion, coagulation, 
and migration show that collisional evolution within a ring or disk of debris leads 
to a few small satellites orbiting Pluto-Charon. These simulations are the first to 
demonstrate migration- induced mergers within a particle disk. The final satellite 
masses correlate with the initial disk mass. More massive disks tend to produce 
fewer satellites. For the current properties of the satellites, our results strongly 
favor initial debris masses of 3 — 10 x 10 19 g and current satellite albedos A pa 
0.4-1. We also predict an ensemble of smaller satellites, R < 1-3 km, and very 
small particles, R ~ 1-100 cm and optical depth r < 10~ 10 . These objects should 
have semimajor axes outside the current orbit of Hydra. 



Subject headings: Planetary systems - Planets and satellites: formation - Planets 
and satellites: physical evolution - Kuiper belt: general 
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INTRODUCTION 



The Pluto-Charon binary is an icy jew el of the solar system. Discovered in 1978 

, and references therein), the central binary 
0.10 ( where M P ess 1.3 x 10 25 g is the mass of 



f lChristv fc Harrington! Il978l : iNoll et al. 



2008 



has a mass ratio, qpc = M c / (Mp + M c ) 



1.5 x 10 24 g is the mass of Char on: iBuie et alJbopd), a n orbital semimajo r 



Pluto and M c 

axis a PC « 17 R P (Pluto radii; 1 R P m 1160 km; lYoung fc Binzellll994!; lYoung et al.l f2007h. 
and an orbital period Tpc ~ 6.4 d (e.g., Sicardy et al. 2006[ Buie et al. 2006 ). Four satellites 
- P5, Nix, P4, and Hydra - o rbit at semimajor axes, qpr ~ 39 Rp, am? ~ 4 3 Rp, ap^ ~ 
49 Rp, and a Hydra « 57 R P flWeaver et al.ll2006l : IShowalter et all 1201 ll . 120121 ). Curiously, 



the orbital periods of the satellites are roughly 3 (P5), 4 (Nix), 5 (P4), and 6 (Hydra) times 
the orbital period of Ch aron. Although t he satellite masses are uncertain by at least an 



order of magnitude (e.g., IBuie et al.ll2006l ; iTholen et al.ll2008l ). they are much smaller than 



the mass of the Pluto-Charon binary. For an albedo A and a mass density p = 1 g cm 
the masses are M P5 «9x 10 16 A" 3 / 2 g, M N ~ ° ' IAl,s 1 ' ~ ° '• ' 1A|7 

and M H « 1,4 x 10 19 A~ 3 / 2 g (e.g., 



2012 



Stansberry et al 



Showalter et al. 2012). For A « 0.04-1 



2008 



2006; 


rti 

Showalter et al. 


2011 




-« ^-.i 
Youdin et al. 


r 


Marcialis et al. 


|l992 


; Roush et al. 


1996|; 



is Af, « 3-300 x 10 19 



Brucker et al.ll2009l ). the combined mass of the four small satellites 
w 2.3-230 x 10" 6 M P . 



tests for the dynamical stability of multiple systems (e.g., ' 




. 3 ealel 


2006; 


T 


lolen et al. 


2008; 


Siili & Zsigmond 


2009; 


Winter et al. 


2010; 


Peale et al. 


2011; 


Youdin et al. 


2012 


). Dv- 



namical fits of orbits to observations yield constraints on the masses, radii, and albedos of the 
satellites. Analyses using the restricted three body problem or n-body simulations identify 
stable orbits surrounding Pluto-Charon and place additional constraints on the masses of the 
satellites. Taken together, current theoretical results for Nix, P4, and Hydra suggest masses 
(albedos) close to the lower (up per) limits derived f rom dynamical fits to observations, with 
M s < 2.5 x 10 20 g and A > 0.2 jYoudin et alJboij ). 



The system also challenges planet formation theories. Currently, most ideas center on a 
giant impact scenario, where the Pluto-Charon binary forms during a collision b etween two 

1 — q)q~ 1 Mp, where q 



1989; 


Canup 


2005. 


2011) 



M P and M 2 



0.25-0.33 (IMcKinnon 



of 5-10 Rp. Over 1-10 Myr, tidal evolution synchronizes the rotat ion of Pluto-Charon 



and then circularizes and expands the binary to its present separation (IFarinella et al.lll979 
Dobrovolskis et al.lll997l ). 



I n some scenar ios, the small satellites also form dur ing the giant impact (e .g. JStern et al. 



20061 ; ICanupll201ll ). Prior to the discovery of P4/P5, IWard fc Canupl (120061 ) developed an 
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elegant model where the satellites move outward during the expansion of the Pluto-Charon 
orbit. If the Pluto-Charon orbit is eccentric throughout the expansion, reson ant interactions 



with C haron can maintain the 1:4:6 period ratio of Charon, Nix, and Hydra. IWard fc Canup 



( 120061 ) comment that a lower mass satellite might orbit with a 5:1 period ratio. However, 
they also note that most of the impact debris probably developed high eccentricity orbits. 
High velocity collisions of debris particles complicate the ability of Charon to maintain a 
3:4:5:6 period ratio for four satellites. 



Lithwick fc Wul (120081 ) consider another path for satellite formation. Their analytic and 
numerical calculations suggest that maintaining the 4:1 (Nix) and 6:1 (Hydra) period ratios 
in an expanding binary requires mutually exclusive eccentricity evolution for the Pluto- 



Charon binary (see also iPeale et al.l 120111 ). Tidal interactions between Nix and Charon 
might establish the 4:1 period ratio, while interactions between Nix and Hydra establish 
the 4:6 period ratio. However, they consider this possibility unlikely. Instead, they propose 
that the Pluto-Charon binary captured many small planetesimals into a disk. Collisional 
evolution then led to the formation of P5, Nix, P4, and Hydra; orbital migration later drove 
the orbits into their current configuration. 

Here, we consider whether the Pluto-Charon small satellites grow approximately in situ 
from material ejected during the impact or captured from the solar debris disk. We begin by 
demonstrating that a giant impact capable of forming Pluto-Charon is a natural outcome 
of scenarios for planet formation at 15-40 AU. After the impact, collisions among smaller 
objects orbiting Pluto-Charon produce a ring of debris. On short time scales, viscosity 
within the ring broadens it into a disk with a radius ad ~ 50-100 Rp. For typical debris disk 
masses Md 10 19 -10 21 g, collisions produce several satellites with masses comparable to 
those of P5-Hydra at semimajor axes reasonably close to the current orbits of the satellites. 
Once formed, satellites may migrate through the leftover debris. These results support the 
idea that satellites with masses comparable to those of P5, Nix, P4, and Hydra grow out of 
material in a circumbinary disk. 

We develop the basic physical processes and calculate relevant time scales for this pic- 
ture in §2. In §3, we describe time-dependent evolutionary calculations of a ring of debris 
surrounding Pluto or Pluto-Charon and its interaction with material orbiting the Sun (§3.1), 
summarize results of coagulation calculations for satellite growth (§3.2), and illustrate plau- 
sible migrat ion scenarios (§3.3). We discuss the implications for the New Horizons mission 



(jSternl 120081 ) and outline how additional calculations can improve and test our formation 
model in §4. We conclude with a brief summary in §5. To isolate basic conclusions from 
analytic estimates and numerical simulations, we use short summary paragraphs and sub- 
sections throughout the text. Table [1] lists the main symbols defined in the text. Table 
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[2] summarizes the physical processes involved in the model and includes a reference to the 
appropriate section of the text. 



PHYSICAL MODEL 



To construct a predictive theory for the formation and evolution of the Pluto-Charon 
binary, we focus mainly on the giant impact scenario. In this picture, two massive protoplan- 
ets within the protoplanetarjj^] debris disk collide and produce a binary planet embedded in 
a ring of debris. Aside from the basic idea of a giant impact, our discussion diverges from 
Earth-Moon models. Unlike the Earth-Moon event, (i) the impact of icy protoplanets does 
not form a disk of vapor and molten material, (ii) an object with roughly the mass of Charon 
survives the impact , and (iii) debris from the impact lies well outside the Roche limit of Pluto 
(jCanupl 120051 . l201ll . and references therein) . 



Our approach has some features in common with recent scenarios where coagulation in 
material just outside Saturn's ring syste m leads to the formation of several of the innermost 



moons (e.g. 



Charnoz et al. 2010 



2011) 



Unlike Saturn, the evolution of the Pluto-Charon 
binary has a large impact on the formatio n and dynamical evolution of the satellites (e.g., 



Ward fc Canupll2006l : iLithwick fc Wull2008l ). The large velocity dispersion of debris from the 
Pluto-Charon impact also complicates coagulation calculations compared to Saturn's rings, 
where the low veloc i ty dispersion allows r apid growth once material escapes the Roche limit 
JSalmon et al.lboiol ; Icharnoz et aDboilj ). 



Our calculations follow some aspects of the IWard fc Canupl (120101 ) method for deriving 
the form ation and evolution of satellites within the circumplanetary di sks of gas giant planets 
fseealso lCanup fc Wardlbood : ISasaki et al.lboid : lOgihara fc Idall2012h . Like lWard fc Canup 
(120101 ). we consider the evolution of a viscous circ umplanetary dis k wher e small particles grow 
into satellites. For satellites around a gas giant, IWard fc Canupl ( 120101 ) derive the evolution 
of a circumplanetary disk of gas and dust fed from the surrounding protoplanetary disk. 
They include the growth and migration of satellites within the gaseous disk. For Pluto- 
Charon, the mass of the gaseous component of the disk is negligible. Thus, we consider 
the evolution of a pure particle disk. Instead of deriving the radial disk structure and the 
growth of satellites analytically, we employ a diffusion code to derive disk structure and we 
use a hybrid coagulation/n-body code to follow the growth and migration of sat ellites as a 
function of radial distance from the central binary (see also lOgihara fc Ida! 120121 ) . 



throughout this paper, we use 'circumplanetary' ('circumbinary') to refer to material orbiting a single 
(binary) planet and 'protoplanetary' to refer to material orbiting the Sun. 
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In this section, we develop the framework for analyzing the formation of the Pluto- 
Charon satellites within debris captured from the protoplanetary disk or retained after the 
giant impact. For the capture hypothesis, we derive estimates for the rate of capture and 
the equilibrium disk mass as a function of the properties of the protoplanetary disk and the 
capture efficiency Our results suggest that capture is a viable mechanism to produce low 
mass satellites surrounding the Pluto-Charon binary. 

For the giant impact hypothesis, we outline the evolution of debris following the impact. 
Using results for orbital stability in binary systems and simple estimates of evolutionary time 
scales for the Pluto-Charon binary and the debris, we conclude that Charon accretes 25% to 
50% of the debris. The binary ejects the rest into a ring surrounding the binary. Because the 
number of 10-100 km objects in the debris is small, a few of them might avoid accretion by 
Charon and end up in the ring. Smaller objects probably suffer a few destructive collisions. 
Thus, the ring consists of, at most, a few large objects and many objects with R < 1 km. 

The evolution of material in the circumbinary ring depends on the relative importance 
of collisional damping, dynamical friction, and viscous spreading. As long as the central 
binary has a n on-zero orbital eccentricity epg, material in the ring has a forced eccentricity 



epc (e.g., 



Heppenheimerlll974j ; iHolman fc Wiegertlll999l ; IWyatt et al.lll999t iLee fc Peak 



20061 ). In a narrow ring, collisions are destructive; the viscous time scale is shorter than the 
time scales for collisional damping and dynamical friction. Thus, the ring spreads into 
a disk. Once the disk has a radial width comparable to its semimajor axis, collisional 
damping becomes more effective. Collisions then produce larger objects. Although viscous 
spreading slows down, the disk continues to expand and may interact with material from 
the protoplanetary debris disk. 

As the largest objects grow, they may migrate through the circumbinary disk. During 
this evolution, gravitational stirring dominates collisional damping. Stirring leads to a colli- 
sional cascade where the smallest objects are ground into small dust grains which are ejected 
from the system. Thus, growing satellites migrate through a disk with a slowly declining 
mass. 

To derive the time scales and important physical issues for this scenario, we begin with 
estimates for the frequency of giant impacts in the protoplanetary disk. After conclud- 
ing that these impacts are common, we demonstrate that a low-mass debris disk probably 
surrounded proto-Pluto prior to impact. We then consider (i) the evolution of debris im- 
mediately after the impact, (ii) the evolution of the expanding debris disk, (iii) interactions 
with the protoplanetary debris disk, and (iv) satellite migration through the circumbinary 
debris disk. 
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2.1. Impact Frequency 

Throughout the history of the solar system, giant impacts are common. In the standard 



mals. small solid obiects with radii of 0.1-100 km (e.e.jGreenbere: et al.lll984t 1 


vVetherill & Stewart 


1989: 


Kokubo & Ida 


19961: IWeidensc 


lilline; et al. 


1997: 


Kenvon & Luul 


1998 


Kenvon 


2002 




Naeasawa et al. 


2005 


Kenvon et al. 


20081). As lone as the mass in planetesimals is larger 



than the mass in protoplanets, dyna mical friction between t he two sets of objects dominates 
viscous stirring among protoplanets (IGoldreich et al.l 120041 ) . Thus, protoplanet orbital ec- 
centricities remain small. Collisions among protoplanets are very rare. When protoplanets 
contain more than half of the total mass, viscous stirring dominates dynamical friction. Pro- 
toplanet eccentricities increase; chaotic growth begins. D uring chaotic growth, protoplanets 



grow through giant impacts with other protoplanets (e.g., lChambersll200ll ; iKominami &: Ida 
20021 : iKenyon &: Bromley! 120061 ) . Giant impacts remain common until protoplanets clear 
their orbits of other protoplanets and leftover planetesimals. 

Outcomes of giant co llisions among protoplanets depend on the impact parame ter, b, and 
the impact velocity Vi mp (lAsphaug et al.ll2006l ; iLeinhardt et al.ll2010t ICanupll201ll ). Usually, 
Vimp > Vesc, where v esc is the mutual escape velocity of a merged pair of protoplanets. Head- 
on impacts with b ?a then produce mergers with little or no debris. Grazing impacts 
with b ~ 1 yield two icy planets on separate orbits and some debris. When b ~ 0.7-0.9, 
collisions often produce a binary planet - similar to the Earth-Moon or Pluto-Charon system 
- embedded in a disk or ring of debris. 

The frequency of collisions that produce a binary planet depends on the distribution of 
possible impact parameters. During the late stages of evolution, the system of protoplanets is 



highly flattened (Wetheril 



19981 : IKenyon fc Bromley 



fc Stewartl 



2006, 



19931 : IWeidenschilling et al.lll997l : IChambers fc Wetheril! 



20101 ) . With most collisions in the plane of the disk, the 
frequency of various types of impacts is roughly the distribution of impact parameters. As- 
suming all impact parameters are equally likely, the probability that a collision will produce 
a binary planet is roughly the probability of an impact with b fa 0.7-0.9, ~ 20%. 



2.1.1. Basic Model 



To estimate the frequency of these impacts among icy objects during the early history of 
the solar system, we examine results from recent numerical calculat i ons of plane t formation 
at 15-150 AU around a solar-type star (IKenyon fc Bromley! 120081 . |2010| . 120121 ) . In these 



calculations, we follow the evolution of protoplanets and planetesimals as a function of 
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semimajor axis a and time t. Calculations begin with an initial ensemble of planetesimals 
with maximum radius Rq and roughly circular orbits in a disk with a surface density of 
solids, 

S(a) = Y, x m a~ n , (1) 

where S is a normalization factor d esigned to yield a sur f ace density eq uivalent to the 
'minimum mass solar nebula' at 1 AU (I Weidenschillingi 1 19 77bl ; lHayashil 1 1 98 ll ) . scaling 
factor to allow a range of initial disk masses, and n = 1-2. 

To derive the outcomes of collisions among planetesimals and protoplanets, we use 
an energy-scaling approach where the outcome depends on the ratio Q c /Q*d, where Q c is 
the collision energy and Q* D is the energy required to disperse half the mass of a pair of 
colliding pla netesimals to infinity. From detailed n-bqdy and smooth particle hydrodyn amics 
calculations (IBenz Asphaudll999l ; iLeinhardt et al.ll2008l ; iLeinhardt k, Stewartll2009i ). 



Q* D = Q b R h + Q g pB? B 



(2) 



where QbR^ h is the bulk component of the binding energy, QgpR 133 is the gravity component 
of the binding energy, and R is the radius of a planetesimal. Table [3] lists the fragmentation 
parameters for the weak, strong, and very strong planetesimals with p = 1.5 g cm -3 used in 
our calculations. 



2.1.2. Results 



Fig. [T] shows the frequency (number per Myr) of impacts with b = 0-1 as a function of 
time for calculations with strong planetesimals (Rq = 1 km) and with several initial surface 
density distributions. In each panel, the curves plot the frequency of impacts between a target 
with mass M t > 10 25 g and a projectile with a mass M p > q p (M t + M p ) where q p = 0.25. 
As summarized in the caption, color indicates the semimajor axis of a disk annulus, ranging 
from 30-36 AU (violet) to 120-150 AU (black). Left-hand panels show results for disks with 
x m = 0.10-0.33 and n = —1.0; right-hand panels plot results for disks with x m = 1.00-3.00 
and n = —1.5. Disk s with (n,x m ) = (—1.0,0.3 3) and (—1.5, 3.00) have the roughly the same 
total disk mass (see Kenyon &: Bromley 201ol ). 



These results demon strate that giant impacts are comm on throughout the evolution of 
the solar system (see also lMcKinnonlll989l ; ICanupll2005l . 120111 ). Initially, all planetesimals are 
smaller than the target mass; the frequency of giant impacts is zero. As planetesimals grow 
into large protoplanets, the frequency of giant i mpacts increases. Because pr otoplanets grow 
faster in the inner disk than in the outer disk ( Kenyon fc Bromley! l2004a| ]bl) . giant impacts 
occur earlier in the inner disk (violet and blue curves) than in the outer disk (magenta and 
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black curves; see also iKenyon fc Bromleyl 120051 ). In each disk annulus, the giant impact 
frequency rapidly rises to a clear peak when the surface density of targets and projectiles 
is largest. As large objects collide and merge into larger protoplanets, the giant impact 
frequency declines. The r ate of decline, p ?: oc t" 1 , is t he rate expected for collisional evolut ion 
in a protoplanetary disk (IKenyon fc Bromleyl 120021 ; iDominik fc Decinl 120031 ; IWyattl 120081 ). 



Several factors contribute to the strong variation of Pi with initial disk mass. In these 
calculations , the number of Pluto-mass obj ects Np grows with increasing disk mass, with 
Np oc x m (IKenyon &: Bromleyl 12008k |2012| ). For giant impacts requiring two Pluto- mass 
objects, the number of possible pairs of impactors scales with the square of the disk mass. 
In more massive disks, the larger mass of leftover planetesimals circularizes the orbits of 
the largest planets more effectively. Thus, gravitational focusing factors also scale with disk 
mass. Combining these two factors, the impact probability grows rapidly with disk mass, 



Pi oc x^, with n 



3-4. 



Finally, the frequency of giant impacts depends on the initial size and strength of plan- 
etesimals (Fig. |2|). In the left panels of Fig. 121 calculations with weak planetesimals yield a 
much lower frequency of giant impacts than calculations with strong planetesimals. When 
planetesimals are strong or very strong , high velocity collisions produce very little debris 
( IKenyon &: Bromleyl 12004a . I201CH . 120121 ). Thus, the mass reservoir to produce Pluto- mass 



protoplanets declines slowly with time. When plane tesimals are weak, collisions produce 



Pi 
2( 



copious amounts of debris ( IKenyon &: Bromleyl 120101 ) ; the mass reservoir declines rapidly 
with time. With less mass to produce massive protoplanets, calculations with weak plan- 
etesimals produce fewer massive protoplanets. Fewer massive protoplanets yield a smaller 
giant impact frequency. 

Although increasing R , the initial maximum size of planetesimals, has little impact 
on the maximum frequency of giant impacts, calculations with larger R produce giant 
impacts later in the evolution than calculations with smaller Rq (Fig. (2J right panels). Large 
planetesimals stir their surroundings more effectively and require long er time scales to grow 
to 1000 km than small planetesimals ( IKenyon fc Bromleyl |2010| . 120121 ) . 



Despite the large range in the giant impact frequency for different initial collisions, all 
calculations yield pi > 0.01 Myr" 1 at 15-30 AU and p { > 0.001 Myr" 1 at 40-60 AU. If the 
distribution of impact parameters is random, the frequency of binary planets is pb p ~ 0.2pj. 
Thus, the frequency of giant impacts capable of producing a binary planet exceeds 1 per 
100-300 Myr at 30 AU. At 20 AU (40 AU), the typical frequency is roughly a factor of two 
larger (smaller). 
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2.1.3. Summary: Giant Impacts are Common 

We conclude that giant impacts capable of producing the Pluto-Charon binary can occur 
at any semimajor axis a ~ 20-40 AU prior to 0.5-1 Gyr. Although calculations starting with 
1-10 km planetesimals produce binary planets earlier in time than calculations with 100 km 
planetesimals, Pluto-Charon binaries are a likely outcome of planet formation in massive 
disks. 

For the rest of this section, we assume a protoplanetary disk similar to the MMSN, with 
n = 3/2 and S = 30 g cm~ 2 . As long as the disk is massive enough to produce a reasonably 
high frequency of giant impacts, other choices lead to similar conclusions. 



2.2. Circumplanetary Debris Disk Evolution: Before and After Impact 



In addition to the two protoplanets involved in the Pluto-Charon- forming giant impact, 
the protoplanetary disk contains a large ensemble of other protoplanets, leftover planetesi- 
mals, and small dust grains. Although protoplanets sweep up this debris inefficiently, mate- 
rial can in teract within the Hill sphere of the protoplan et and form a circumplanetary debris 



disk (e.g.. IWeidenschillingl |2002| ; iGoldreich et al.ll2002l ). Once a circumplanetary disk forms 



mate rial orbiting the planet may interact with other material passing through the Hill sphere 



(.e.g. 



Durda fc Sterol l2000l : ISteroll2009l : IPires dos Santos et al.ll2012f ). 



The Hill radius, Rh, defines a spherical surface where the gravity of a planet roughly 
balances the gravity of the Sun. For a planet with mass density p orbiting the Sun at 
semimajor axis a : 



i? H ft;4x 10 e 



R 
Rp 



P 



2 g cm 



1/3 



a© 



20 AU 



km . 



(3) 



When an object from the protoplanetary disk enters the Hill sphere, there are four possible 
outcomes: 



It can interact with Plut o-Charon, find a temporary orbit, and collide with other 



objects on similar orbits ( IPires dos Santos et al. 



20121 ). Because the collision time is 



much longer than the lifetime in a captured orbit, captured material probably cannot 
form a circumplanetary disk. Thus, we do not consider this option. 



It can collide with another object from the protoplanetary disk (e.g.. IWeidenschilling 
20021 ; iKominami et al.ll201ll ). This collision yields a single large object and some debris. 
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If this material can lose enough kinetic energy to achieve an orbit with semimajor 
axis a < ^Rh and 7 ~ 0.4, it may become bound to the protoplanet. Solids with 
7 > 0.4 cannot execute closed orbits around the planet and are eventually ejected (e.g. 
Martin fc Lubowlhoilf ). 



It can collide with an object in the circumplanetary disk. During the late stages of 
planet formation, incoming projectiles have much larger velocities than target objects 
orbiting within the disk. If the projectile has a mass much smaller than the target, 
the target - a nd other nearby objects within t he disk - accrete most of t he incoming 
material (e.g.. iDurisen et al.lll989t ISternll2009l ; iHousen fc Holsappldl201ll ). Projectiles 
with masses much larger than the target remove material from the disk. 

It can pass through the Hill sphere without any collision. When the optical depths of 
the circumplanetary and protoplanetary disks are smaller than unity, this outcome is 
the most likely one. 



For simplicity, we first assume that material entering the Hill sphere at a distance 
a < 7-R# from the planet can interact with other objects and remain bound. This material 
encounters the Hill sphere at a rate Md pa Z'kYiqVLq^Rh) 2 '', where S is the surface density 
of the protoplanetary disk, Q Q is the angular velocity of the planet around t he Sun, and the 
factor of three accounts for a s mall degree of gravitational focusing^ (see also Weidenschilling 
20021 ; iKennedy fc Wyattll201ll ). Adopting conditions appropriate for a proto-Pluto, this rate 
is 



M^5x lO 21 / x 



''0 



g yr 



(4) 



,20 AIL 

where we include a factor f Q to account for time-dependent depletion of the initial surface 
density distribution. 

The probability that an object in the Hill sphere interacts with other objects is roughly 
Td + t p , the sum of the optical depths of the circumplanetary and the protoplanetary disk. 
With r pa E/-R, we estimate 



t p «3x 10" 



i f<3 x r 



Rr, 



100 km 



3/2 

20 AIL' 



(5) 



where R ma x is the radius of a typical large object in the protoplanetary disk. 



For the circumplanetary disk, the optical depth for particles moving through the Hill 
sphere is pa (£/ R)(H z /a), where H z = v z VL~~ l is the vertical scale height and v z is 



2 During this epoch, the vertical scale height of the disk, Hq, is much larger than Rh ■ 
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the vertical velocity dispersion. Usually, v 7 is ha lf the radial velocity dispersion v r (e.g., 
Hornung et al.lll985l ; IOhtsukilll992l ; llda et al.lll993[ ) . The radi al velocity dispersion is roughly 



the escape velocity of the largest solid objects in 1 


he disk 


1995: 


Kenvon & Bromlev 


2001: 


Goldreich et al. 


2004). 


and mass density p ps 1 g cm 3 , 



(e.g., 



Barge fc PeUatlll990l : lKokubo & Ida 



53 



R„ 



cm s 



(6) 



1 km / 

The optical depth in the circumplanetary disk for collisions with material in the protoplan- 
etary debris disk is then independent of particle size: 

M d 



r d ^ 10 



-10 



10 20 gj V20 AU 



(7) 



where M d is the total mass of the disk. 



To estimate the probability that collisions within the Hill sphere yield some material 
bou nd to the planet, we first consider the fate of two objects from the protoplanetary disk (see 
also IWeidenschilring||2002l ; iGoldreich et al.ll2002l ). Most material interacts near the outer edge 
of the disk at ^Rh where the orbital velocity around the planet is v or b ~ 2000 cm s _1 . With 
typical impact velocities exceeding ~ 2 x 10 4 cm s _1 , a single object must lose roughly 99% of 
its ki netic energy to r e main bound. In his analysis of the formation of binaries in the Kuiper 
belt, IWeidenschillingl (120021 ) derives capture probab ilities of 0.1% to 1%. For conditions in 
a young protoplanetary disk, IGoldreich et al.l (120021 ) explore whether dynamical friction or 
three-body interactions allow capture. They conclude that ~ 3% of incoming material can 
remain bound to the protoplanet. At times when giant impacts are common, dynamical 
friction is more than a factor of ten less effective. Reducing the impact of dynamic al friction 
implies a capture probability of < 0.3%, comparable to the IWeidenschillingl ( 120021 ) estimate. 



For incoming planetesimals colliding with material already in the circumplanetary disk, 
depletion is much more likely than capture. Because the disk is mainly collision debris, 
incoming planetesimals are much larger than objects within the disk. Typical collisions 
between an incoming planetesimal and an orbiting target then produce an impact crater on 
the planetesimal and some ejecta from the target. For an incoming planetesimal with R m 
10-100 km and a target wit h 0.1-1 km, the mass in th e ejecta is less than a few per cent 
the mass of the target (e.g., iHousen fc Holsapplell201ll. and references therein). Althoug h 



circumplanetary objects may accrete the ejecta (e.g., iSternl l2009i : IPoppe fc Horanyill201ll ). 
the incoming planetesimal escapes with most of the mass of the target and depletes the 
circumplanetary disk. 



When inco ming objects have small masses , collisions add material to the circumplane- 
tary disk (e.g.. iDurisen et al.lll989l : ISternl 120091 : IPoppe fc Horanyill201ll ). However, objects 
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with R < 0.1 km contain much less than 1% of the mass of the protoplanetary disk. Thus, 
incoming low mass objects add relatively little mass to the circumplanetary disk. 

In steady-state, depletion from high velocity collisions with the disk balances the ad- 
dition of debris from collisions within the Hill sphere. Defining as the fraction of debris 
captured by the planet, steady-state implies fdT p t a . Solving for the state-state disk mass: 

fd \ ( Rmax \ 1 / «0 \ 3 / 2 /0 n 



M d „ 3 x 10-/^ (Jl) (J^J ( 



20 AU 



Viscosity within the disk probably transports mass inward onto the central planet and out- 
ward through the Hill sphere. Thus, this disk mass is a rough upper limit. With a typical 
input rate of roughly 10 11 g yr" 1 , it takes at least 30-100 Myr to reach this steady-state disk 
mass. 

As long as the protoplanetary disk has a surface density reasonably close to its ini- 
tial value, Pluto-mass planets can sustain modest circumplanetary disks of solid material. 
This mass is roughly at the lower limit for the sum of the masses of the known satellites 
surrounding Pluto-Charon Growth of pla nets or satellites is rarely 100% efficient (e.g., 



Kenyon fc Bromley! l2004al . 12008k |2010| . 120121 ) . Unless capture is more efficient than our sim- 
ple estimates, forming massive satellites around Pluto-Charon requires an additional source 
of material. The giant impact which produces Charon is the most likely source of this 
material. 



2.3. Immediate Aftermath 



After the impact, the Pluto-Charon binary and the surrounding debris evolve on various 
time sca les. Initially, the binary has semimajor axis ape ~ 5-10 Rp and eccentricity epc ~ 
0.1-0.8 (jCanupll201ll ). The initial orbital period sets the shortest time scale in the system: 

3/2 



PC 



ape 
5 R P 



d. 



(9) 



Debris from the impact mostly lies in an elliptical ring along Charon's orbit. For a 
binary wi th mass ratio q = 0.1, th e minimum semimajor axis for stable orbits outside the 
binary is ( IHolman fc Wiegertlll999l ) 



1.96 + 4.68e - 2.17e 2 PC )a PC . 



(10) 



For epc 
inside a n 



0.1-0.8, test particle orbits with a min < 2.4-4.3 a PC are unstable. Any debris 
is accreted by the binary pair or ejected into orbits outside the binary. Material 
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outside the binary likely collects in a ring with orbital semimajor axis > a m i n . This 
material has a minimum orbital period 

\ 3/2 

'.Amin « 4- 9 1^1 d. (11 ; 



5 Rj 

If the debris from the impact is distributed fairly smoothly along the orbit, Charon 
accretes this material at a rate Mq ~ T,flpcRc ~ 0.1 Ma Tpc~ l , where Rc is the radius 
of Charon and flpc is the angular velocity of the Pluto-Charon binary. The time scale to 
sweep up all of the debris is 

t acc = M d /M c w 10 J d . (12) 

The fraction of debris particles accreted by Charon depends on the time scale to eject 
material from the vicinity of the binary. Dynamical simulations of particles along Charon's 
orbit suggest it takes a few orbital periods to accrete or to eject particles from the vicinity of 
the binary system. Thus, Charon probably accretes 25% to 50% of the debris along its orbit. 
The Pluto-Charon binary ejects the rest. To conserve energy and angular momentum, the 
Pluto-Charon binary contracts. For a mass in debris, M d ~ 10 20 g at a min , the decrease in 
the Pluto-Charon orbital period is much smaller than 1%. 

As Pluto-Charon eject particles beyond a min , collisions among debris particles help to 
circularize their orbits. If the debris maintains a fairly smooth surface density along the 
orbit, the collision rate is dn/dt pa (Md/7r5a PC )Qpc / pR, where 5 is the width of the ring 
and R is the radius of an object in the debris. With p = 1 g cm -3 , the collision time is 

tc = w 3 (>) («) (^) (PS. V' 2 d . (13 ) 



x o.iy vo.ikmy \ M d J \5 R_ 

For a mass, M d pa 10 20 g, the collision time for 10 km and larger particles exceeds the time 
scale for Charon to eject or to sweep up the debris. Thus, large proto-satellites can escape 
the Pluto-Charon binary more or less intact. 

While Charon ejects or sweeps up the debris, smaller particles suffer destructive colli- 
sions. Around the Pluto-Charon binary, typical orbital velocities are 0.5 km s _1 . With initial 
e = 0.1-0.8, impact velocities are roughly 100 m s -1 . Impact energies exceed 10 7 erg g -1 . Col- 
lisions involving particles with R ~ 0.3-5 km are completely disrupted; among 10-30 km ob- 
jects, cratering co llisions reduce the mass by roughly 5% (30 km) t o 25% (10 km) per collision 



see, for example, iLeinhardt et al.ll2010l ; lKenyon fc Bromleyil2012l ). For M d pa 10 20 g and the 
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collision time in eq. f lT3|) . mutual collisions among particles with R < 5 km (M <5x 10 17 g 
for p « 1 g cm -3 ) produce an ensemble of smaller particles orbiting the Pluto-Charon binary. 

As Charon ejects debris, collisional damping effectively circularizes the orbits of the 
smallest objects. For objects with R < 0.1 km, the damping time is compara ble to or smaller 



than the orbital period of the Pluto-Charon binary (e.g., eq. [131 see also iGoldreich et al. 



20041 1 . Efficient damping leaves the debris on roughly circular orbits with semimajor axes, 



Because collisions among large objects are infrequent (eq. [13]), their orbits are circu- 
larized by dynamical friction or by direct erosive collisions with the smaller objects (e.g., 



Goldreich et al.ll2004f ). Following the impact, particle velocities exceed the escape velocity 
of 10-30 km objects. Thus, dynamical friction is in the dispersion regime where gravita- 
tional focusing factors are roughly 1. The damping time scale for dynamical friction is then 
tdf ~ (R/H a )Tpc, where S s is the surface density of the small particles. If small particles 
contain a fraction f s of the initial mass, the damping time scale is 

5 \ ( R \ /10 20 g\ ( a PC ^ 7/2 



*"Mi»Aan=J bsrJ m) d (14) 

Dynamical friction circularizes the orbits of large particles on long time scales compared to 
the orbital time of the binary. Because collisions with small particles erode the large particles 
on similar time scales, large particles may not survive the circularization process. 

Interactions with any pre-existing circumplanetary debris disk also circularize the orbits 
of impact debris. Dynamical friction and collisional damping reduce the radial velocity 
dispersion of the impact debris. With a typical mass of ~ 10 19 g spread over 1000 Rp, the 
pre-existing disk has an angular momentum roughly 10 times larger than the initial angular 
momentum of the impact debris. Circularizing the orbits of impact debris beyond a m j n 
requires a modest reduction in the angular momentum of pre-existing disk material. The 
time scale to circularize the orbits of the 0.1-1 km particles in debris is roughly the collision 
time, ~ 0.1-1 yr. 

We conclude that the Pluto-Charon binary ejects debris into a ring with a typical 
semimajor axis > a min . For simplicity, we set ~ 4apc ~ 20Rp. The solid material 
in the ring consists primarily of small particles with R < 1 km and perhaps a few larger 
objects with R > 10 km. Th e evolution of the central binary then depends on tidal processes 



between Pluto-Charon (e.g.. lDobrovolskis et al.lll997l ) and gravitational interactions between 



Pluto-Charon and the surrounding ring of debris. 
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2.4. Orbital Evolution of the Central Binary 

Over long time sca les, tidal interactions bet ween Pluto and Charon m odify the orbit and 
rotation of the binary ( IDobrovolskis et al.lll997l ). Following iPeald (119991 ) . the rate of change 
of the angular rotational velocity, u, of Charon is to ~ (— 45/76)(pf2 4 i? 2 7 //iQ)sign(u; — Q), 
where 1/Q is the angular phase lag in the tidal oscillation and fi is the rigidity. Holding the 
structure of Pluto-Charon fixed at Q ~ 100 and fi ~ 4 x 10 10 dyne cm -2 , the time scale to 
synchronize u is 

'^y^-^„, us) 

s a PCfi J \t r ,o t r J 

where t r is the rotational period, and t r ^ is the initial rotational period. Thus, Charon's 
rotation synchronizes rapidly with the orbit on a time scale, t sync ~ 10 2 -10 4 yr for apc~ 
5-10 R P . 

Tidal interactions also expand the orbit. Tides change the angular orbital frequ ency at 
a rate ft « (-9/2)(kM c Rpn 16/3 )/ (OM P \G( M P + Mn)} 5 / 3 ) signfo;-^ JPeale! ll999h. where 
k w 9.4 x 10~ 3 is the Love number for the potential ( Love 1944 ; Peale 19991 ). If the structural 
parameters are held fixed, tides slowly expand the Pluto-Charon semimajor axis: 



t PS 

h sync ~ 



170 



ape — apc,o 1 



texp J 



6/39 



(16) 



where the expansion time is t exp 6500(apc,o/5-Rp) 39//6 yr. ^From an initial semimajor axis 
apcp ~ 5-10 Rp, the orbit evolves to the c urrent a PC = 17 Rp in 0.2-20 Myr (see also 



Farinella et al.l Il979l ; IDobrovolskis et al.l 119971 ) . 



Dynamical interactions with the circumbinary debris disk help to circularize the Pluto- 
Charon orbit, but they slow orbital expansion. For a ring of material with ~ 20Rp 
surrounding a Pluto-Charon binary with ape ~ 5Rp and epc ~ 0.1, numerical simulations 
with the Orchestra code (see §3.3) suggest 



and 



(hi 



-2.5 x 10" 



-2 x 10" 



Md 



10 20 g 



M, 



yr 



(17) 



R P yr" 



'11 



10 20 g ^ 

When the debris disk has a mass comparable to the masses of the satellites, the time scale 
for dynamical interactions to change a and e is similar to the tidal time scale. 
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2.5. Initial evolution of debris ring 



2.5.1. The ring spreads 

Once material lies in a ring outside the Pluto-Charon binary, collisional processes within 
the ring are more important tha n the dynamical evolution of the bina ry. Tidal forces from 



the binary truncate the disk (e.g.. lLin fc Papaloizoulll979l; |Pringldll99lf) and impose a forced 
eccentricity, e/ ~ epg, on ring particles (e.g., Murray fc Dermottl 1999 : Wyatt et al. 1999 



Lee &: Peald 120061 ; iTsukamoto fc Makinol 120071 ). The time scales for these interactions (eqs. 
[T3] and [14]) are much smaller th an the 10 5 -10 7 yr time scale for tidal interactions to 
circularize and to expand the orbit (IDobrovolskis et al.lll997l ). 



Viscous spread ing occurs on a time scale t v ~ (5a) 2 / v where v is the viscosity (ILynden-Bell fc Pringle 



1974 ; |Pringldll98ll ). For a disk of identical indestructible pa rticles undergoing inelastic col- 
lisions, the viscosity is roughly ( IGoldreich fc Tremaindll982l ) 

v 2 0.46r 

(19) 



where r 



n i + t 2 

3E/ 4Rp is the optical depth of the ring. 



When r is small, gray itational scattering augments the collisional viscosity (IShu &: Stewart 



1985 



1985|) 



Hornung et al 



Il985h . 



The viscosity from gravitational encounters is roughly (IHornung et al. 



0.5 



v 2 v 2 



n 



V 



(20) 



r 



where v esc is the escape velocity of the particles and L = _ha[H z /2R) 3-5. In an equilibrium 
disk, v r < v esc ( jHornung et al.l Il985l ; iBarge fc Pellatl Il990h . The gravitational scattering 
com ponent of the viscosi ty then exceeds the collisional component by a factor of 3-10 (see 
also 



Hornung et al.lll985l ). 



Initially, particles in the ring have large velocity dispersions set by the dynamics of the 
central binary. The gravitational component of the viscosity is then much smaller than the 
collisional component. With v r ~ epc Qa^, the viscous time during the early parts of the 
evolution is 

R \ ( M d \ f_e_\ ( <i,i 
0.1, 



t, 



650 



0.1 km/ ViO 20 g7 V0.1/V 2 0^p 
With an initial viscous time of 1-2 yr, the ring spreads rapidly 

The ratio of the collision time to the viscous time is 



7/2 



(21) 



Q 5a<i 



(22) 
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For a narrow ring with v r ~ epc Qa>d and badjdd ~ epc*, the collision and viscous times are 
roughly equal. Each particle suffer several destructive collisions as the ring spreads. 



2.5.2. Particle sizes in the ring 

When destructive collisions dominate, small particles are ground down into smaller and 
smaller objects. Radiation pressure can eject very small particles with R w 1-100 /mi on 
the local orbital time scale; Poynting- Robert son drag pulls particles into the central star or 
planet on much longer time scales. For particle orbits around the S un, a balance betw een 



radiation pressure and gravity establishes an absolute minimum size (IBurns et al.lll979l ): 



0.6 jum 

where p is the mass density of a particle and (3 is the ratio of radiation pressure to gravity. 
Solar radiation pressure ejects particles with (3 > 0.5-1; thus R m i n ~ 0.5-1 pm for icy 
particles with p w 1 g cm -3 . 

For orbits around a planet, the Sun e jects particles when radiation pressure exceeds 



the gravity of the planet (IBurns et al.lll979l ). For Pluto-Charon, the minimum particle size 



increases by more than an order of magnitud^E 



1/2 '20 AU X 1/2 



iU,«7«20 pm, (24) 

\AapcJ \ a J 

where a is the orbital semimajor axis of Pluto-Charon around the Sun. For the Pluto- 
Charon binary, particles with R < 0.5-1 pm are ejected from the binary and the solar system. 
Part icles with R 1~20 urn are ejected from the binary but remain bound to the Sun (see 



also iPoppe fc Hor anvil l201lh 



Despite the increased efficiency of radiation pressure in removing small particles from a 
circumplanetary disk, Poynting- Robertson drag is ineffective. For particles with p rs 1 g cm -3 
on Sun-centered orbits, the time scale for Poynting-Robertson drag to pull particles radially 
inward is 

^ ( R \ ( a & 



•"""""i.i^teJ"' (25) 



3 We assume radiation pressure ejects particles when /3 ~ v/vq, where v is the orbital velocity around the 
planet and Vq is the orbital velocity of the planet around the Sun. Numerical simulations suggest this limit 
is appropriate for a disk of small particles where collisional damping counteracts radiation pressure. 
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For orbits around Pluto-Charon, the drag time is roughly 30% longer (IBurns et al.lll979l ). 
Small particles with stable orbits around Pluto-Charon, R > 20 /im, spiral toward the 
central binary on time scales comparable to the time scale for tidal evolution of the central 
binary. These particles probably cannot pass inside the tidal boundary of the binary at 
roughly a > 2apc- Thus, small particles are likely swept up by large planetesimals close to 
the central binary. 

Initially, the largest object in the ring has a radius R < 10-30 km (m < 3-100 xlO 18 g). 
Disruptive collisions can reduce these objects in mass by a factor of ten. Thus, the maximum 
size is roughly 1-10 km. If all large objects suffer destructive collisions, the maximum size is 
close to 0.1 km. Near the binary, the smallest particles have radii of roughly 20 /im. In this 
size range, 20 fim < R < 0.1 — 10 km, collisional evolution probably produces a differential 
size distribution n(r) oc r -3 5 (e.g.. lDohnanyilll969l ; IWilliams &: Wetherilllll994j . but see also 



Belyaev & Rafikov 2011). Thus, most of the mass is in the largest objects. 



2.5.3. After the ring spreads 



As the ring spreads, collisional damping becomes more effective. The radial velocity 
dispersion declines. In this regime, the Hill radius, Rh, sets the length scale. Scaling eq. ([3]) 
to the size of the Pluto-Charon binary, 



R 



H 



10 



R 



1 km 



20 Rp 



km 



(26) 



If a few large particles with R > 10 km survive dynamical ejection into the ring, they have 
Hill radii of 0.1-0.3 Rp. These objects st i r up material within a few Hill radii, ~ 0.3 — 1 Rp, 



of their orbits (e.g., Ilda fc Makinolll993l : lOhtsuki et al.l 120021 ) . Because the initial width of 



the ring is a few Pluto radii, a few large particles can stir the entire ring. Although small 
particles with R < 0.1 — 1 km have much smaller Hill radii, they are very num erous. In a ring 
with mass pa 1Q 20 g, small p articles can also stir up the entire ring (e.g., iKokubo fc Ida 



19951 : iKenvon fe Bromlevl 120011 ). 



When particle stirring dominates collisional damping, the minim um radial velocity dis 



1990 



persion is roughly the escape velocity of the largest objects, eq. (El) ( iBarge fc Pellat 
Kokubo fc Idalll995l ; IKenvon fc Bromley! l200ll ; iGoldreich et al.l I2004J ) . Small ring particles 



have vertical velocity dispersions of roughly half the radial velocity dispersion, v z ~ v r /2. 
The vertical scale height of the ring, H z = v z Q _1 , is: 



30 



Rmax 

1 km 



20 R F 



3/2 



km . 



(27) 
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The typical vertical height of the particle disk is smaller than a Pluto radius. 

When Scid/cid w 1 and v r ~ tv,mm> the ratio of the collision time to the viscous time is: 

If no large fragments with 72 > 10 km survive the giant impact, growth is faster than viscous 
spreading. The typical growth time is 

R \ /I0 2 °g\ I a d 



( < K 50 litoj \mf) y- ^ 

Thus, the largest particles grow rapidly to 10-30 km. 

After the largest objects reach radii of 10-30 km, the time scales for growth and viscous 
spreading are roughly equal. As viscosity continues to expand the ring into a disk, the largest 
objects continue to accrete smaller debris particles. Eventually, the disk consists of a few 
large objects orbiting the Pluto-Charon binary. 



2.6. Migration Through the Circumbinary Disk 



Although viscous spreading and collisional processes probably dominate the early evo- 
lution, gravitational scattering plays an important role in the local structure of the disk. 
Large objects try to clear their orbits by scattering smaller objects. Viscous spreading 
among smaller objects tries to fill the orbits of the larger objects. When the large objects 
dominate, they reduce the surface density along their orbits and increase the surface density 
outside their orbits. Because the surface density enhancements are not smooth, the large 
objects feel a torque from the small objects. If the sum of all of the torques does not vanish, 
the large objects migrat e radially inward or outward (e.g.. lGoldreich fc Tremaindll982uWard 
19971 : iKirsh et al-lbopgh . 



When scattering overcomes viscous spreading, large objects can open gaps in the radial 
distribution of planetesimals. For convenience, we define the Hill radius necessary for a 
particle to op en up a gap i n a dynamically cold disk . For the Pluto-Charon binary, this 
radius is (e.g., lRafikovll2001t iBromley fc Kenyonl 120131 ) 

1/3 / p \ -1/3 / N 1/3 

lo ( ^— r ) | I (-?-) [^-] km (30) 



1 L gap 



53 cm s 1 ) ^lO 20 g / 
15 km have physical radii R 



1 km 



20 R P 

26]). For circumbinary disks 



Objects with Rh ~ 15 km have physical radii R ~ 1 km (eq 
with Md pd 10 20 g, satellites with R > 1 km open gaps in the disk. All four Pluto-Charon 
satellites have R > 3 km. Thus, each can open a gap in a cold circumbinary disk. 
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Fo r sate llites with R > R gap , there are two possible modes of migration (e.g.. IBromley fc Kenyon 



201 lbl . 120131 ) . If the satellite can clear a gap in its corotation zone and migrate across this 
gap in one synodic period, the satellite undergoes fast migration. Satellites with Hill radii 
in the range R gap < Rh < Rfast satisfy this condition. For the Pluto-Charon binary, 



R 



fast 



20 



10 



20 



1/2 



20 Ri 



3/2 



km. 



(31) 



Depending on the disk properties and satellite location, this Hill radius corresponds to objects 
with physical radii of a few kilometers. Coupled with the limits on R gap , satellites with 
physical radii of 1-5 km undergo fast migration. These objects drift radially inward or 
outward at a rate 

a f^t ~ ±5 ( tt^- 



3/2 



km yr 



-1 



(32) 



10 20 gj \20Rp 

With minimum radii of 12 km (Nix) and 15 km (Hydra), the two massive satellites of Pluto- 
Charon are too massive to undergo fast migration. If the smaller P4 and P5 have albedo 
A > 0.5, they are small enough to migrate in the fast mode. 

Massive satellites that open up a gap in the disk migrate through the disk at a rate 
that depends on the disk viscosity. In gaseous disks with a large viscosity, this type II 
migration can transport gas giants from 5-10 AU to within a few stellar radii of a solar - 



type star in ~ 1 Myr (ILin fc Papaloizoul Il986bl ; IWardl 119971 : iNelson fc Papaloizou 



In a disk of particles, the smaller viscosity leads to a smaller type II rate (llda et al. 



Cionco fc Bruninil |2002| ; iKirsh et al.l 120091 ; iBromley fe Kenyon 



disk in Pluto-Charon, the expected migration rate is (IBromley fc Kenyon! 120131 ) 



2004 



2000; 



2011bh. F or a c ircumbinary 



4 „ „ -0.15 (-%-) («) (^V /2 km yr- . 



(33) 



Satellites migrate more rapidly (i) when the disk is more massive, (ii) when the satellite is 
more massive, and (iii) when the satellite is closer to the central binary. For standard disk 
and satellite masses, the Pluto-Charon satellites migrate roughly 1 Rp every 10 3 yr. 

The ability of satellites to migrate efficiently depends on the ratio of the vertical scale 
height H z to the Hill radius Rh- When collisional damping dominates particle stirring, 
H z /Rh ^ 1- Satellites can open gaps in the disk and migrate at the nominal rates. In 
a disk dominated by particle stirring, H Z /R H & 3(a/20R P ) 1 ^ 2 (see eqs. [2T)H2"T] ). When 
H z /Rh > 1, satellites cannot open gaps in the disk. Migration is then 'att enuated,' with a 
rate smaller by a factor of roughly (HjR H f w 30-100 at a w 20-60 R P llda et alJlioOO 



Kirsh et al.l 120091 ; IBromley &: Kenyon! 1201 lbl ). Despite this reduction, the attenuated type 
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II rate is still substantial, ~ 1 Rp every 10 Myr, similar to the time scale for tides and 
dynamical interactions to modify the semimajor axis of the binary orbit. 

This discussion suggests two plausible phases for migration in a circumplanetary disk 
surrounding Pluto-Charon. During the early stages of the evolution, rapid migration of small 
satellites is possible when collisional damping dominates particle stirring. Later on, when 
stirring by newly-formed satellites dominates collisional damping, attenuated migration can 
modify satellite orbits as the inner binary expands. 



2.7. Summary of the Physical Model 

Our estimates suggest that satellite growth and migration are likely within a circumbi- 
nary debris disk produced from material (i) accreted from the protoplanetary disk or (ii) 
surviving the Pluto-Charon giant impact. If the disk is massive enough, satellite formation 
in captured material is straightforward: collisional damping produces a disk where small 
particles can grow into satellites. In debris from a giant impact, material surrounding a 
Pluto-Charon binary with an initial semimajor axis ape = 5Rp (see also Table [2]) follows 
the following evolutionary scenario. 



The i mpact produces the Pluto-Charon binary with debris along Charon's orbit (jCanup 



201ll ). Charon either accretes or ejects the debris on a time scale t ~ 1-3 d. 



Particles that survive accretion/ejection lie within a circumbinary ring with semimajor 
axis ad ~ 20Rp. As they make their way into the ring, all debris particles suffer several 
destructive collisions. Because they are so rare, larger particles might escape the binary 
more or less intact. 

Shortly after the ring forms, the time scale for viscous spreading is roughly 10 yr. 
Thus, the ring spreads into a disk. 

Once the disk has 8a /a ~ 1, collisional damping promotes the growth of larger objects 
within the spreading disk. Material from the protoplanetary debris disk adds material 
to the disk, enhancing the growth rate. 

When particles have sizes of a few km, they migrate through the leftover smaller 
particles. Although longer than the growth rate, the migration rate is sufficient to 
transport large objects throughout the circumbinary disk in less than 0.1 Myr. 
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NUMERICAL CALCULATIONS OF SATELLITE FORMATION IN THE 

PARTICLE DISK 



Although we cannot perform a complete numerical calculation of satellite formation 
following a giant impact, we can divide the problem into several pieces which illustrate the 
likely outcomes. After considering the viscous evolution of a circumbinary disk of particles, 
we examine the growth and the migration of satellites within this disk (see also Table |2]) . 

To perform these calculations, we use Orchestra, an ensemble of computer codes for 
the formation and evolution of planetary systems. Currently, Orchestra consists of a ra- 
dial diffusion code which computes the time evolution of a gaseous or particulate disk 
(IBromley &: Kenyonll2011al ). a multi annulus coagulation code which c omputes the time evo- 
lution of a swarm of planetesimals (IKenyon fc Bromley! l2004al . 120081 ). and an n-body code 
which comp utes the orbits of gravitationally interacting protoplanets (IBromley &: Kenyon 
20061 . l2011al ). Within the coagulation code, Orchestra in cludes algorithms for treating inter- 



action s between small particles and the gaseous disk (e.g.. lAdachi et al 



1976: 



Wei denschilling 



1977a|) and between coagulation mass bins and n-bodies (jBromley fc Kenyonll2006l ). To treat 
interactions between small particles and the n-bodies more rigorously, the n-body calcula- 
tions include tracer particles. Massless tracer particles allow us to calculate the evolution 
of the orbits of small particles in response to the motions of massive protoplanets. Massive 
tracer particles enable calcu lations of the response of n -b odies to the changing gravitational 
potential of small particles ( Bromley fc Kenyon 2011al lbl. 12013 ). 



3.1. Viscous Evolution 

We begin our set of numerical calculations with the evolution of rings of solid particles 
orbiting a proto-Pluto or the Pluto-Charon binary. The ring has initial mass Mj, semimajor 
axis and width Sad- To explore a broad range of plausible initial conditions, we first 
consider the evolution of a ring of particles orbiting a central object with the mass of Pluto. 
Next, we examine the evolution around a central mass equivalent to the Pluto-Charon binary. 
Although not precisely equivalent to a ring orbiting an eccentric binary, these calculations 
show how different initial conditions change the global evolution of the ring. 



-23 - 



3.1.1. Diffusion Calculations 



To follow the evolution of th e surface density of a ring su r rounding a p l anet, we solve 



the radial diffusion equation (e.g., iLynden-Bell &: Pringldll974l ; |Pringldll98ll . Il99lf ) 



— = 3a- 1 — 
dt da 



,1/2 







i~ {^« V2 } ) - ( ^ 



dE\ 



(34) 



where a is the radial distance from the central star and t is the time. The first term represents 
the change in surface density E from viscous evolution; the second term allows changes in E 
from external processes, including material added from the protoplanetary debris disk. 



To derive the appropriate viscosity for the disk, we follow ISalmon et al.l (120 10f ) and 
ignore the gravitational scattering component of the viscosity. Although s cattering dom- 
inate s when disk particles are in equilibrium (e.g., eqs. (19-20); see also IShu &: Stewart 
19851 ; iHornung et ajJ[l985J ), collisions dominate the viscosity during the early stages of the 
evolution. Thus, our approach yields a good estimate for the viscous time early on, at the 
expense of a factor of 2-3 underestimate at late times. For this initial exploration of disk 
evolution in Pluto-Charon, this accuracy is acceptable. 

When the disk is gravit at tonally stable, there are two sources of viscosity from collisions: 
random motions of particles within the disk {y t ) and sound waves traveling between colliding 
particles (z/ c ). When the self-gravity of the disk exceeds the rotational shear and the internal 
velocity dispersion, self-gravity adds another component (v g ). Self-gravity is important when 
Q T < 1-2, where 



Qi 



Oil 



(35) 

and G is the gravitational constant ( Toomrel 1964c Goldreich &: Lynden-Bell 1965 ). The 
collisional viscosity is independent of Qt ( Araki fc Tremainel 19861 ): 



R 2 Qt 



(36) 



For icy particles with p = 1 g cm -3 , r = 3E/4i2. Thus, v c grows linearly with the particle 
size and the surface density. The viscosity drops with increasing distance from the central 
object. The transl ational viscosity depends on the Toomre stability parameter Q T (e.g., 
Daisaka et aJboOlh : 

G 2 S 2 



13^ 

o,2 



'h 
.46t^ 



Qt < 2 
Qt >2 



(0ASt\ 
K 1+r 2 ) 

where is the ratio of the Hill radius to the radius of a pair of identical particles, 
Rh/2R- The self-gravity component of the viscosity is 

v t Q T < 2 

Q T > 2 



(37) 
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(38) 
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From §2.2, the mass input to the disk is 



— = 3E fi (/r p - T d ) , 



(39) 



ext. 



where / ~ 10 3 (§2.2); r p (eq. [5]) and r d (eq. [7]) are the optical depths of the protoplanetary 
and circumplanetary disks. 

To solve eq. (|34|) . we substitute x = 2a 1//2 and S = xE. On a grid with 1001 grid 
points uniformly spaced in x from the surface of Pl uto to OARh, we use a standard explici t 



method to derive the evolution of S with time (e.g. JBath fc Pringldll98ll : iPress et al.lll992l ). 
Calculations with a fully implicit approach yield identical results. For explicit and implicit 
approaches, the calculations conserve mass and angular momentum exactly. For disks sur- 
rounding single and binary planets, we adopt a flow through outer boundary condition, with 
E = at a = OARh. For a single planet, we adopt the flow through boundary condition 
at the planet' s surfa ce, a = Rp. To treat truncation of the disk by the inner binary, we 
follow iPringld ( )199ll ) and set the radial velocity of disk material at this boundary to zero. 
Formally, the boundary condition is 



0. 



(40) 



We apply this boundary condition at a 



(eq. [TO 



3.1.2. Before Impact: Circumplanetary Disk Evolution 

To explore whether circumplanetary disks have sufficient mass to produce the satellites 
orbiting Pluto-Charon, we consider the evolution of a low mass ring surrounding a Pluto- 
mass protoplanet. As a standard model, we adopt small particle sizes, 10 3 cm, for material 
in the circumplanetary disk, which implies a small radial velocity dispersion, 0.5 cm s _1 . 
The calculation begins with a low mass (Ma = 10 11 g) ring at a = 20Rp. To derive the 
maximum disk mass, material is added at the rate from eq. ( 139]) with a depletion factor / 
= 1. 

Fig. [3] illustrates the long-term growth of the circumplanetary disk. New material 
from the protoplanetary disk gradually erases the initial conditions and produces a roughly 
constant surface density from the surface of the proto-planet to the edge of the stable region 
at 0.4 Rh- As the system evolves, the surface density slowly increases. After an evolution 
time of 100 Myr, the system begins to reach a steady-state with a total disk mass of Md^s ~ 
2.5 x 10 19 g. Although the steady-state surface density distribution peaks close to the central 
proto-planet, the distribution is shallow, with E oc a~ n and n w 0.20-0.25. 
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The behavior of the surface density in Fig. |3]has two main causes. Early on, the disk has 
little mass and very low optical depth. Material from the protoplanetary disk rarely interacts 
with the circumplanetary disk. Thus, there is little collisional depletion. During this period, 
the surface density remains roughly constant with radius. As the evolution proceeds, the 
disk optical depth grows, eventually enabling collisions with disk material. Because the disk 
has a larger scale height at the outer edge of the disk, there is relatively more depletion 
from the outer disk than the inner disk. This extra depletion produces a maximum surface 
density close to the central planet. 

Throughout the evolution, the disk transfers material inwards onto the central planet 
and outwards into the protoplanetary disk. After ~ 1 Gyr, the planet accretes a negligible 
amount of material, ~ 10 12 g. Because most disk material is added near the outer radius, 
the disk transports a much larger amount of material into the protoplanetary disk, ~ 10 14 g. 
This material makes a negligible contribution to the ~ 10 25 g orbiting within 5-10 Hill radii 
of the planet. Thus, material from the circumplanetary disk has a limited impact on the 
mass of the central planet or the evolution of the protoplanetary debris disk. 

The mass of the steady-state circumplanetary disk scales with the adopted accretion 
efficiency /, the depletion of the protoplanetary disk / Q , and the radius of a typical large 
object in the protoplanetary disk R ma x- The numerical results are close to the analytic 
estimate in §2.2: 

M d „, » 2.5 x !0>% (^) (fe) (^y 3 ' 2 • (41) 

Although the steady-state disk mass is independent of the particle size in the disk, 
the amount of material accreted by the central planet and the amount transported out 
through the disk vary with particle size. In the outer disk, larger viscosity enables more 
mass transport. Thus, mass loss through the outer edge of the disk grows roughly linearly 
with the particle size. In the inner disk, mass accretion depends on the viscosity and on the 
rate angular momentum flows out through the disk. Larger viscosity increases the efficiency 
of both processes. Thus, the mass accretion rate grows somewhat less rapidly with particle 
size. In all cases, the total mass accreted onto the planet or lost into the protoplanetary disk 
is smaller than the steady-state mass in the disk. 

3.1.3. After Impact: Circumbinary Disk Evolution 

We now consider the viscous evolution of a circumbinary ring of material formed during 
the Pluto-Charon giant impact. The calculations assume a narrow ring of solid particles with 
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total mass Md at initial semimajor axis a = 20 Rp. The particles have radius R = 0.1 km 
and radial velocity dispersions v r = 5, 50, or 500 cm s _1 . In the first set of calculations, the 
Pluto-Charon binary has fixed separation ape = 5 Rp. In the second set of calculations, 
the binary expands from an initial ape = 5 Rp to the current ape = 17 Rp over several Myr 
(eq. [IE]). 

Fig. H] shows results for calculations with = 3 x 10 20 g. Initially, the solids lie in 
a narrow annulus at a = 20Rp, roughly half the semimajor axis of the current satellites. 
At the start of the calculation, the ring spreads roughly uniformly inward and outward. 
Eventually, the inner edge of the disk reaches a semimajor axis a m i n = lapc- The inner 
boundary condition prevents material from flowing in past this semimajor axis. The inner 
edge then remains fixed while the outer radius spreads outward. 

At all times, the extent of the disk depends on the radius and the radial velocity dis- 
persion of the solids (Fig. H]). In disks with R = 0.1 km and v r = 5 cm s _1 (blue curves), 
viscous spreading is slow. It takes roughly 1 Myr for the inner edge of the disk to reach 
a = lapc and more than 10 Myr for the outer edge of the disk to reach the current orbits 
of the small Pluto-Charon satellites, a ~ 40 — 60Rp. Disks with larger velocity dispersions 
evolve more rapidly. When v r = 500 cm s _1 (maroon curves), the inner edge of the disk 
reaches a ~ lapc in only 10 3 yr. At ~ 10 4 yr, the outer edge of the disk crosses the orbits 
of P5, Nix, P4, and Hydra. 

Far away from the inner and outer edges, the surface density follows a power law with 
£ oc a~ n . During the early stages of the evolution, the power law is steep, with n pa 2. As 
the evolution proceeds, the surface density distribution becomes shallower, with n pa 1.0-1.5. 
After roughly 1 Gyr, the outer edge of the disk reaches a pa OARh- The slope is then very 
shallow, with n pa 0.5-1.0. 

When we include the tidal expansion of the binary orbit, the evolution of the solids 
changes dramatically (Fig. [5]). As the binary expands, the inner boundary condition main- 
tains the inner edge of the disk at a = 2apc- Thus, the inner disk expands in step with 
the binary orbit. Throughout the disk, the viscosity is independent of the inner boundary 
condition. The expanding binary has little impact on the evolution of the outer edge of the 
disk. With a larger inner disk radius and a similar outer disk radius, the disk covers a much 
smaller radial extent compared to calculations without an expanding inner binary. 

Despite differences in the radial extent of the disk, the evolution of the slope of the 
surface density distribution is fairly independent of the evolution of the central binary. Early 
on, £ oc a~ n with n pa 2. As the binary approaches its current separation of 17 Rp, the 
slope slowly declines to n pa 1.0-1.5. At very late times, ~ 100 Myr to 1 Gyr, the slope is 
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shallower still, with n ~ 0.5-1.0. 

Calculations with different combinations of disk mass, particle size, and radial velocity 
dispersion lead to similar results. The time scale for the outer edge of the disk to reach the 
current orbits of the satellites depends on the initial mass in the debris, the typical radius of 
particles in the disk, and the radial velocity dispersion of disk particles. For our calculations, 
the time scale for outer disk material to reach a = 50 Rp is 

Massive disks composed of small particles with large velocity dispersions evolve on the fastest 
time scales. 

The scaling in eq. ( pf2l) depends solely on the translational viscosity v t . For M r i fh 
10 18 — 10 22 g, v t contributes more to the viscosity than v c (see eqs. |36]-[37]). Because 
v t oc Sf 2 i? _1 , the viscosity grows linearly with increasing disk mass, falls linearly with 
increasing particle size, and grows with the square of the radial velocity dispersion. Thus, 
the viscous time scale t u oc RM^ x v~ 2 . 

For the satellites orbiting Pluto-Charon, these results have several clear consequences. 



Because the disk and the inner binary expand at comparable rates, understanding the 
origin of the 3:4:5:6 period ratio for the satellites requires treati ng the evolution of the 



central binary and th e disk of debris simultaneously (see also IWard fc Canupl 12006 



Lithwick & Wu! 120081 ). 



In §2, our results suggest that debris ejected from the vicinity of the Pluto-Charon 
binary probably consists of many small particles with R < 0.1 km and perhaps a 
few large objects, R ~ 10-20 km, which do not suffer destructive collisions during 
the ejection process. The diffusion calculations indicate that debris with a few large 
objects - and hence large velocity dispersions - evolves more rapidly than debris with 
no large objects. Thus, understanding the details of satellite formation requires better 
knowledge of the evolution of debris during the impact. 

Disks with smaller radial velocity dispersions are narrower in semimajor axis than disks 
with larger radial velocity dispersions. If migration is unimportant, the current orbits of 
P5, Nix, P4, and Hydra are easier to understand if the initial velocity dispersion of the 
debris is small. Rapid migration allows a broader range of initial velocity dispersions. 
We return to this point in §3.3. 
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• Disks with smaller radial velocity dispersion expand out to the current orbits of P5, 
Nix, P4, and Hydra on much longer time scales than disks with larger radial velocity 
dispersion. If satellites form rapidly, on time scales less than 10 4 — 10 5 yr, the orbits of 
the satellites depend on how they respond to the expansion of the inner binary and the 
expansion of the disk. If the satellites form on 1-10 Myr, interactions with the remnant 
disk are more important than interactions with the central binary. We consider the 
formation time in §3.2. 



3.1.4- Summary of Viscous Disk Calculations 

The results of the diffusion calculations show that the satellites orbiting the Pluto- 
Charon binary are somewhat more likely to form in the debris from the giant impact than 
in material accreted from the protostellar disk. Unless accretion from the protostellar disk 
is more efficient than our estimates suggest, the mass of a circumplanetary disk surrounding 
a proto-Pluto is a factor of 1.5-10 times smaller than the current mass of the satellites. Th e 



amount of debris produced in a giant impact is uncertain by a factor of 10-100 (ICanupll201ll ). 
Some impacts yield masses large enough to produce P5, Nix, P4, and Hydra. 

Although debris from the impact initially lies well inside the current orbits of the satel- 
lites, viscous diffusion and tidal evolution of the Pluto-Charon binary drive disk material 
to larger orbits on interesting time scales. Expansion of the binary prevents disk material 
from spreading inside a 2apc- Once Charon achieves its current orbit with ape = 17 Rp, 
disk material orbits outside a ~ 34Rp. Thus, the dynamics of the binary orbit precludes 
satellites inside the orbit of P5. Viscous spreading transports disk material to large radii, 
a > 60i?p, in roughly 1 Myr. Because a broad range of M d , R, and v r yield disk material at 
the current orbits of the four satellites, the satellites could plausibly form within a variety 
of disks. 

To place additional constraints on satellite production, we now examine the growth of 
satellites in a disk of particles. 



3.2. Growth of Satellites in the Particle Disk 

To calculate the growth and evolution of satellites in a disk of small particles, we use a 
hybrid multiannulus coagulation- iV-body code. We compute the collisional evolution of an 
ensemble of particles in a disk orbiting a central object of mass Mpc- The code uses sta- 
tistical algorithms to evolve the mass and velocity distributions of low mass objects with 
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time and an iV-bo d y alg orithm to follow the individual trajectories of massive objects. 
Kenyon fc Bromlevl (120081 ) describe the statistical (coagulation) code; iBromley &: Kenyon 
(120061 . l2011al ) describe the iV-body code. Here, we briefly su mmarize the basic aspec ts of 
our approach (for a general review of collisional evolution, see lYoudin &: Kenyonll2012l ). 



We perform calculations on a cylindrical grid with inner radius a« n and outer radius 
a out . The model grid contains N concentric annuli with widths 5cii = 0.025aj centered at 
semimajor axes eij. Calculations begin with a cumulative mass distribution n c (m ik ) oc m~ k q 
of particles with mass density p = 1 g cm -3 and maximum initial mass m . For comparison 
with investigators that quote differential size distributions with n oc r~ q , q' — q + 1. Here, 
we adopt an initial q' Q = 0.17; thus most of the initial mass is in the largest objects. 

Planetesimals have horizontal and vertical velocities h ik (t) and fjfc(t) relative to a 
circular orbit. The horizontal velocity is related to the orbital eccentricity, ef k (t) = 1.6 
{hik{t)/VK,i} 2 , where Vk,% is the circular orbital velocity in annulus i. The orbital inclination 
depends on the vertical velocity, if k {t) = sm~ 2 (2(vik(t)/VK,i))- 

The mass and velocity distributions of the planetesimals evolve in t ime due to inelastic 
collisions, dra g forces, and gravitational encounters. As summarized in Kenyon &: Bromley 
(j2004al 120081 ) . we solve a coupled set of coagulation equations which treats the outcomes of 
mutual collisions between all particles with mass rrij in annuli aj. We adopt the particle- in- a- 
box algorithm, where the physical collision rate is navf g , n is the number density of objects, 
a is the geometric cross-section, v is the relative velocity, and f g is the gravitational focusing 
factor. Depe nding on physical conditio ns in the disk, we derive f g in the dispersion or the 
shear regime (IKenyon fc Bromley! 120121 ). For a specific mass bin, the solutions include terms 
for (i) loss of mass from mergers with other objects and (ii) gain of mass from collisional 
debris and mergers of smaller objects. 

Collision outcomes depend on the ratio Q c /Q*d, where Q* D is the collision energy needed 
to eject half the mass of a pair of colliding planetesimals to infinity and Q c is the center 
of mass collision energy. To derive Q* D (eq. [2]), we adopt parameters for weak particles 
(Table [3]). For two colliding planetesimals with masses mi and m 2 , the mass of the merged 
planetesimal is 

m = mi + ^2 — m>ej j (43) 
where the mass of debris ejected in a collision is 



m 



0.5 (m 1 + m 2 ) (j^j 



9/8 



(44) 



To place the debris in our grid of mass bins, we set the mass of the largest collision frag- 
ment as itll = 0.2m esc and adopt a cumulative mass distribution n c oc m~ qd with q^ = 
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0.833, roughly the value expected for a system in collisional equilibrium ( Dohnanvi 



1969 



Williams fe Wetherilllll9~9~4l : kanaka fc Idalll996l : b'Brien fc GreenbereJl2003l ; kobavashi fc Tanaka 



2010l ). This approach allows us to derive ejected masses for c atastrophic collisions with 



Qr ~ Q*n and for cratering collisions wit 



Williams fc Wetherill[[l99i kanaka fc Ida! 



i Q r < Q* n (see also IWetherill fc Stewart! 



O'Brien fc Greenberg 



2003; 



1996 



Stern fc Colwelllll997l: iKenyon fc Luu 



Kobavashi fc Tanaka 



19931 : 
19991 : 



2010|). 



To compute the evolution of the velocity distribution, we include collisional damping 
from inelastic collisions and gravitational interacti ons. For i nelast ic an d elastic collisions, w e 
follow the statistical, Fokker-Planck approaches of lOhtsukil (Il992l ) and lOhtsuki et al.l ( 120021 ). 
which treat pairwise inter actions (e.g., dynam i cal fri ction and viscous stirring) between all 
objects in all annuli. As in lKenyon &: Bromley! ( 120011 ) . we add te rms to treat the probability 
that objects in annulus %\ interact with objects in annulus i<± (IKenyon fc Bromlevi l2004al . 



2008 ). We also compute long-range stirring from distant oligarchs ( Weidenschilling 19891 ) 



In the iV-body code, we directly integrate the orbits of objects with masses larger than 
a pre-set 'promotion mass' m pro . The calculations allow for mergers among the iV-bodies. 
Additional algorithms treat mass accretion from the coagulation grid and mutual gravita- 
tional stirring of n-bodies and mass batches in the coagulation grid. To treat interactions 
among proto-satellites, we set m pro = 10 17 — 10 19 " 



Our calculations follow the time evolution of the mass and velocity distributions of 



objects with a range of radii, r, 
ejected by radiation pressure, we set r m 
than the largest object in each annulus. 



r min to T; t j 



r max- To track the amount of material 
20 ^m. The upper limit r max is always larger 



3.2.1. Pure Coagulation Calculations 

To explore possible outcomes for an ensemble of solid particles orbiting the Pluto- 
Charon binary, we consider a suite of pure coagulation calculations. The calculations begin 
with a ring of material surrounding a central mass equivalent to the Pluto-Charon binary. 
The total mass in the ring ranges between Md = 3 x 10 19 g and Md = 3 x 10 21 g. The ring 
consists of 32 concentric, circular annuli extending from 15 Rp to 50 Rp. This cylindrical 
grid provides a reasonably accurate representation of orbits of debris particles from the giant 
impact around a circular binary system with an initial separation of 5 Rp. 



Debris captured from the protoplanetary disk lies at much larger disk radii then con- 
sidered in these calculations. Because coagulatio n calculations are easily scaled to different 
semimajor axes ( IKenyon fc Bromlevi 120081 . 120101 ). it is straightforward to apply our results 
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to the larger disks produced by capture. 

Within each annulus, the solid particles range in size from 20 /mi to 0.1 km. Consistent 
with solutions for the radial diffusion equation in §3.1, the initial surface density distribution 
is £ oc a -1 . To allow the system to find an equilibrium between the large orbital eccentricity 
of material ejected from the central binary and the small eccentricity expected from collisional 
damping, we adopt an initial eccentricity eo = 0.1 and an initial inclination io/ e o — 0.5. 

All of these calculations follow a similar pattern. With the large initial (e, i) of all 
particles, collisions among large particles are destructive and produce copious amounts of 
smaller particles. Among the smaller particles, collisional damping is effective. Through 
dynamical friction, smaller particles damp the orbits of the larger particles. Thus, the 
(e, i) for all particles declines with time. Collisions become less destructive. Eventually, 
collisional velocities become small enough to promote growth over destruction. During this 
phase, roughly 25% of the initial mass is ground down into particles smaller than 20 /im. 
This material is ejected from the binary. 

After the system reaches a low velocity equilibrium, the largest particles begin to grow 
rapidly. During a very short phase of runaway growth, objects grow from roughly 0.1 km to 
3-30 km. As the largest objects grow, they gravitationally stir the smaller particles. Once 
the largest objects contain most of the mass, stirring dominates collisional damping. The 
orbital (e, i) of the smaller objects rapidly increases, reducing gravitational focusing factors 
and halting runaway growth. At the end of runaway growth, each annulus contains a handful 
of proto-satellites with radii of 3-30 km. 

Following runaway growth, the evolution of the system stalls. The proto-satellites 
rapidly accrete any small particles left in their annuli. Because the small particles con- 
tain so little mass, this growth of the large particles is very small. In the particle-in-a-box 
approximation adopted for these pure coagulation calculations, the probability of collisions 
among large objects is small. Thus, the sizes of the proto-satellites reach a rough steady 
state set approximately by the initial mass in each annulus of the grid. 

Figure M illustrates the growth of satellites for calculations with = 3 x 10 19 g (thin 
lines) and Md = 3 x 10 21 g (thick lines). As summarized in the caption, each line represents 
the evolution in the radius of the largest object in annuli ranging in distance from 15 Rp 
to 50Rp from the center-of-mass. In these disks, satellite growth is very rapid. Despite 
an initial period of destructive collisions, it takes only 1-100 yr for an ensemble of 0.1 km 
fragments to grow into 10 km satellites. The growth time is inversely proportional to the 
initial disk mass, t g w 100(3 x 10 20 g/Md) yr. Once satellites reach a maximum mass of 
5-50 km, growth stalls. Satellites remain at a roughly constant radius for thousands of 
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years. 

These calculations demonstrate the inevitable growth of 5-30 km satellites from an en- 
semble of much smaller objects. Independent of the initial orbital eccentricity, collisional 
damping reduces the internal velocity dispersion to levels that promote mergers during col- 
lisions. Despite the loss of roughly 25% of the initial mass to a collisional cascade, small 
objects grow rapidly into satellites with radii of 5-30 km. The growth time is short com- 
pared to the time scale for viscous spreading of the disk (Fig. [5]) or the time scale for binary 
expansion (eq. [IS]). Thus, the growth of small satellites orbiting the Pluto-Charon binary 
is very efficient. 

Despite the efficiency of small satellite formation, coagulation calculations produce too 
many satellites. In these models, each of the 32 annuli in the calculation yields 0-3 satellites 
with radii of 5-30 km. Even with small number statistics, this result is much larger than 
the current number of known satellites. However, the coagulation calculations do not allow 
large-scale dynamical interactions among satellites. For the models shown in Fig. [6j satellites 
in adjacent annuli have orbital separations of 5-15 Rp. Thus, satellites are close enough to 
perturb the orbits of their near est neighbors, leading to ch aotic interactions as in simulations 



of terrestrial planet formation (IKenyon fc Bromley! 120061 ) . 



3.2.2. Hybrid Calculations 

To investigate dynamical interactions among newly-formed satellites, we now consider 
a suite of hybrid simulations with Orchestra. Objects in the coagulation code which reach 
a preset promoted into the n-body code. When a satellite with (m^-, e^, z^) 

in an annulus with semimajor axis a; is promoted, we assign (e n , i n ) = (e^iy), a n = aj + 
gAcii, and a random orbital phase, where g is a random number in the range —0.5 to +0.5. 
The n-body code converts these orbital elements into an initial (x, y, z) position and an 
initial (v x ,v y ,v z ) velocity vector. The n-body code follows the trajectories of all promoted 
satellites. Algorithms within Orc hestra allow n-bodies to interact with coagulation particles 



(Bromley fc Kenvonll2006l . l2011ah . 



The starting conditions for these calculations are identical to those for the pure coagu- 
lation models. We consider three cases for the initial disk mass, = 3 x 10 19 g (low mass 
disk), ma = 3 x 10 20 g (intermediate mass disk), and ma = 3 x 10 21 g (massive disk). For 
each case, we adopt a different promotion mass: m pro = 3 x 10 17 g (low mass), m pro = 10 18 g 
(intermediate mass), and 10 18 g (high mass). These promotion masses are a 

compromise between the 'optimal' promotion mass, m pro = 3 — 30 x 10 16 g, required to allow 
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newly formed n-bodies to adjust their positions and velocities to existing n-bodies and the 
practical limits required to complete calculations in a reasonable amount of time. For each 
combination of disk mass and promotion mass, we ran 20-25 simulations. 

To check the accuracy of our results for the intermediate and high mass disks, we ran 
an additional 12 simulations for each of two cases with promotion masses a factor of three 
larger and a factor of three smaller than the promotion masses listed above. When the 
promotion mass is a factor of three larger, the chaotic growth phase begins later and lasts 
longer. Often, chaotic growth in the outer disk is (unphysically) well-separated in time from 
chaotic growth in the inner disk. This feature inhibits scattering and radial mixing among 
proto-satellites. When the promotion mass is a factor of three smaller, the onset, character, 
and duration of chaotic growth changes very little. In all cases, the final number of satellites 
is fairly independent of promotion mass. 

Because fragmentation has a small impact on the results of pure coagulation calcula- 
tions per unit cpu time, these hybrid calculations do not include fragmentation. Neglecting 
fragmentation increases the mass available for satellite growth by roughly 25%. Without 
fragmentation, collisions fail to produce copious amounts of 1 m to 100 m particles. Colli- 
sional damping among these debris particles and dynamical friction between the debris and 
larger survivors reduces orbital e and i, aiding runaway growth. Thus, these hybrid calcula- 
tions artificially increase the evolution time required for satellites to reach R > 10 km. 

Fig. [7] shows the evolution of the semimajor axes for promoted satellites orbiting a single 
central object with the mass of the Pluto-Charon binary. All objects begin with r < 0.1 km 
(m < 4 x 10 12 g). The evolution time required for these objects to reach the promotion mass 
scales with the initial disk mass, t pro 10 3 (3 x 10 19 g/m^) yr. This time scale is roughly 
a factor of ten longer than derived for pure coagulation calculations with fragmentation. In 
the inner disk, objects have shorter orbital periods and shorter collision times. Thus, the 
first promoted objects appear near the inner edge of the disk. As the calculation proceeds, 
satellites are promoted farther and farther out in the disk. Eventually, a few promoted 
objects appear in the outer disk. 



Every calculation experiences chaotic growth (jGoldreich et al.ll2004l ; iKenyon &: Bromley 



20061 ). During chaotic growth, the satellites scatter one another throughout the grid. Because 
objects grow fastest in the most massive disks, chaotic growth begins earlier - ~ 100 yr - 
in the high mass disk and much later - ~ 10 4 yr - in the low mass disk. Chaotic growth 
is also 'stronger' in more massive disks. In more massive disks, massive large satellites are 
more numerous and generate larger radial excursions of smaller satellites than in lower mass 
disks. 
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Throughout the evolution shown in Fig. [7J the n-bodies slowly accrete the very small 
planetesimals remaining in the coagulation grid. The number of leftovers correlates well with 
the number of n-bodies. Thus, the total mass in leftover planetesimals within the coagulation 
grid approaches zero as the number of n-bodies reaches a minimum. 

In all of these calculations, the final number of satellites correlates inversely with the 
initial disk mass. In massive disks with = 3 x 10 21 g, there are usually 2 or 3 massive 
satellites with r rj 50-80 km at the end of the calculation. As we reduce the initial disk 
mass, the calculations produce more satellites with smaller masses. In intermediate mass 
disks with Md = 3 x 10 20 km, chaotic growth leaves all of the initial mass in 4-5 satellites 
with radii of 20-30 km. In low mass disks, there are very few mergers during chaotic growth. 
After ~ 10 7 yr, there are typically 7-9 satellites with radii of 7-12 km. 

To test how outcomes change with the initial size distribution, we consider a suite of 
15 simulations with two additional 10 km planetesimals placed randomly in the grid. These 
large planetesimals stir their surroundings and thus counteract collisional damping by nearby 
small planetesimals. With much larger masses than any other planetesimals in the grid, they 
have large collisional cross-sections and can grow very rapidly. 

Fig. [8] compares the time evolution of the semimajor axes for n-bodies in calculations 
with (lower panel) and without (upper panel) two additional massive planetesimals. In 
standard calculations with = 3 x 10 19 g, there is a short period of chaotic growth at 
~ 10 4 — 10 5 yr followed by a long quiescent period with an occasional merger. When there are 
two large planetesimals at the onset of the calculation, these planetesimals slowly accumulate 
small planetesimals. Meanwhile, small planetesimals in the rest of the disk rapidly grow into 
sizes (4-5 km) that allow promotion into the n-body code. From Fig. [8], the time scale for 
these objects to reach the promotion mass is roughly 10 3 — 10 4 yr independent of the two 
large satellites. However, stirring by the two large satellites promotes an earlier oligarchic 
growth phase which enables more uniform growth of the largest planetesimals. Thus, there 
are many more promoted objects and a more intense chaotic growth phase. In addition to 
the two initial massive planetesimals, one or two other planetesimals accrete all of the other 
promoted objects. After ~ 10 5 yr, there are only a few small n-bodies left. The 3-4 massive 
satellites accrete these objects in a few hundred thousand years. 

This result is typical of all calculations in low or intermediate mass disks that begin 
with a few large (10 km) planetesimals placed randomly throughout the grid. Stirring by 
the large objects slows runaway growth and allows a larger group of oligarchs to grow rapidly. 
Although these oligarchs reach the promotion mass on similar time scales, many more of them 
reach the promotion mass. Typically, calculations with initial large planetesimals produce 
2-4 times as many n-bodies as calculations without large planetesimals. Stirring by the 
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10 km planetesimals and interactions among the many n-bodies leads to a very active phase 
of chaotic growth. During chaotic growth, the radial excursions of the n-bodies are 2-3 times 
larger in semimajor axis. Collisions among n-bodies are more frequent, leading to a system 
with fewer, but more massive, satellites. With a set of two initially large planetesimals, low 
(intermediate) mass disks produce 3-4 (2-3) massive satellites instead of 7-9 (4-5). 

For all of these hybrid calculations, the radius evolution follows a standard pattern 
(Fig. [H]). As in Fig. [6J the largest objects in the coagulation grid grow from 0.1 km to 
1-10 km. Once they are promoted into the n-body grid, they interact with objects in 
the coagulation grid and all of the promoted objects in the n-body grid. Unlike the pure 
coagulation calculations, growth does not stall at 5-30 km. The extra interactions between 
promoted objects and the larger volume sampled by scattered objects enables a few large 
objects to accumulate all of the remaining mass in the grid. Thus, these objects reach radii 
that are 30% to 100% larger (and 2-8 times more massive) than the largest objects in the 
pure coagulation calculations. 

To conclude this section, Fig. [TU] plots pairs of (e, i) for satellites that survive for 10 7 yr. 
For most satellites, the orbital inclination is small: the average inclination is < % > = 0.1° ± 
0.2°. Despite the apparent excess of high i objects at small e for the low and intermediate 
mass disks, the average i and the dispersion in the average do not vary with initial disk mass. 
The average e grows with disk mass, however, with < e >= 0.02 ± 0.02 for M d = 3 x 10 19 g, 
< e >= 0.03 ± 0.02 for M d = 3 x 10 20 g, and < e >= 0.06 ± 0.06 for M d = 3 x 10 21 g. For 
small and intermediate initial disk masses, the average eccentricity is very low. For large 
disk masses, the < e > is a factor of 2-3 larger. 

Although the final orbital elements of satellites are independent of the promotion mass, 
e depends on the initial masses of the largest planetesimals. When we add several large 
planetesimals to the initial distribution of 0.1 km and smaller planetesimals, the calculations 
yield fewer satellites with larger masses (Fig. [8]). As in calculations with larger total disk 
masses, a few large satellites stir up their surroundings more than many smaller satellites. 
These satellites then have larger e. 

3.2.3. Summary of Coagulation and Hybrid Calculations 

The suite of pure coagulation and hybrid calculations demonstrates that the Pluto- 
Charon satellites grow rapidly from a disk of small planetesimals. When calculations include 
fragmentation, destructive collisions grind down roughly 25% of the initial disk mass into 
20 /im particles which are ejected from the binary system. Collisional damping among 
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larger debris particles reduces orbital eccentricity of 0.1-10 m objects. Dynamical friction 
between these small planetesimals and the surviving large planetesimals damps the orbital 
eccentricities of the largest objects, greatly reducing the frequency of destructive collisions 
and enabling rapid growth of surviving planetesimals. Thus, fragmentation removes mass 
and aids the eventual rapid growth of satellites. 

In hybrid calculations without fragmentation, less collisional damping and dynamical 
friction slow growth considerably. However, satellites still grow on time scales similar to the 
expansion time for the central binary. The number of satellites N s , typical satellite radius 
R s , and the time for satellites to reach their final mass tf scales with the initial mass of 
the planetesimal disk. For calculations without any large (10 km) fragments at t — and 
M d « 3 - 300 x 10 19 g, we infer 

N s « 2 22 - 5 " log Md , (45) 



for the number of satellites 



for the typical radius, and 



for the growth time. 



fi,^27.5[' 3X ^° Mg '| km, (46) 



In low and intermediate mass disks, calculations with several 10 km fragments produce 
fewer satellites with larger masses. For w 3 — 30 x 10 19 g, N s ps 3-5 and R s ~ 15-60 km 
instead of the N s ps 4-10 and R s ps 7-30 km for calculations without large fragments. Be- 
cause large fragments tend to sweep up any small planetesimals, calculations with fragments 
yield a smaller spread in the number and masses of planetesimals than calculations without 
fragments. 

The orbital elements of model satellites agree reasonably well with the observed elements 
of the Pluto-Charon satellites. Roughly 10% of the model satellites have e < e#, where e# 
is the orbital eccentricity of Hydra. Nearly all (83%) satellites have orbital inclinations 
i i$ i-H, where in is the inclination of Hydra. Satellites with integer period ratios are also 
common. For calculations with M d = 3 x 10 21 g, 60% of the outer satellites have period 
ratios reasonably close to the 5:3, 2:1, or 3:1 commensurabilities with the inner satellite. 
Disks with smaller initial masses produce more closely-packed, lower mass satellites which 
often lie close to the 3:2 and 4:3 commensurabilities. For these lower mass disks, 50% to 
70% of satellites have period ratios within a few per cent of the 4:3, 3:2, 5:3, 2:1, or 3:1 
commensurabilities. 
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These results are encouraging. Model satellites often lie close to the observed orbital 
commensurabilities in the Pluto-Charon satellites. The derived inclinations of model satel- 
lites also agree very well with observations. Although lower mass disks produce satellites 
with the smallest e, derived eccentricities for all calculations are factor of 4-20 larger than ob- 
served. Thus, the calculations have two successes - the small inclinations and the likelihood 
of close orbital commensurabilities - and one failure - the lack of very small eccentricities 
for the satellites. 

To explore the outcomes of formation models in more detail, we now consider numerical 
calculations of satellite migration through the circumbinary disk surrounding Pluto-Charon. 
Aside from testing the analytic estimates for R gap , Rfast, and a from §2.5, these calculations 
allow us to estimate the frequency of migration-driven mergers and to examine the evolution 
of e and i for an ensemble of migrating satellites. 



3.3. Migration of Satellites Through the Circumbinary Disk 



To explore satellite migration, we use the iV-body component of the Orchestra code. In 
this mode, the calculations directly track interactions between satellites and massive "super- 
particles" that represent small particles in the disk. For massive obje cts within a disk around 



a single central object, our calculatio ns (Bromley Sz Kenvon 2011a b[) reproduce published 



resul ts of previous investigators (e.g., Malhotra 19931 ; Hahn fc Malhotra 1999 ; Kirsh et al. 



20091 ). Our investigation of migration within Saturn's rings includes extensive tests with gaps 
in the disk, orbital resonances, and the gravit ational perturbations of distant massive moons 
outside resonance (IBromley fc Kenyonl 120131 ). To test the theoretical limits on migration 
rates in a disk surrounding Pluto-Charon, we first consider calculations of a lone satellite 
with R w Rgap — Rfast within a particle disk around a single central object with the combined 
mass of Pluto-Charon. We then consider how a binary central object modifies the mode and 
rate of migration. Because the hybrid calculations often produce many small satellites, we 
conclude this section with simulations of multiple satellites orbiting a central binary. 

For this suite of simulations, we adopt a disk surface density distribution, E oc a -1 , with 
a fixed mass of 3 x 10 20 g. The disk extends from a = 20 Rp to a = 70 Rp around a single 
object with a radius of 1 Rp or a binary with a separation of 5 Rp. We follow an ensemble 
of super-particles and satellites in an annulus of full width of 5 Rp. Super-particles have 
masses of l/2000 th the mass of the satellite. These objects interact with the satellite and the 
central mass, but not with each other. Satellites have fixed bulk mass density, p — 1 g/cm 3 . 



At the start of each simulation, super-particles are dynamically cold. Thus, the disk 
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is geometrically thin, with H z /Rh < 1. Particle trajectories evolve solely by interactions 
with the central (single or binary) object a nd massive satellites. Unlike our simulations of 
Saturn's A ring (IBromley fc Kenyonl 120131 ). there is no collisional damping among super- 
particles. These initial conditions are ideal to assess a satellite's ability to migrate through 
the disk and lead to fairly robust upper limits on the migration rate. In a more realistic disk, 
collisional damping, dynamical friction, and viscous stirring generally produce increases in 
the velocity dispersion and vertical scale height of disk particles on time scales comparable to 
the growth time (see §3.2) . Because migration rates fall off as the inverse cube of the m ean 
disk particle eccentricity (llda et al.ll2000l ; iKirsh et al.ll2009l ; IBromley fc Kenyonl l2011bl ). we 
expect that migration time scales in a real disk are somewhat longer than our estimates for 
idealized disks. 



3.3.1. Migration of a Single Satellite 

Figure [TT1 illustrates results for single satellites around a single central mass (left panels) 
and a central binary system. At a = 20 Rp (Fig. [ITJ upper left panel), satellites with R « 1- 
4 km should open a gap and undergo fast migration (eqs. [30H3T] ) . Satellites with R = 3 km 
migrate at rates a ~ 20 km yr -1 , close to the rate predicted in eq. f )32|) . As the satellite 
radius increases from 3 km to 5-10 km, the migration rate drops considerably. These larger 
satellites have too much inertia to maintain the fast radial drift rate. Larger satellites with 
R > 10 km migrate at rates close to predicted rates for type II migration (eq. [33]). 

Migration rates also depend on the semimajor axis of the satellite (Fig. [TTJ lower left 
panel). When the velocity dispersion of the disk is fixed, it is harder for smaller satellites 
with R = 3-4 km to open up a gap in the disk at larger a than at smaller a (eq. [3U] ). Thus, 
small satellites migrate much more slowly at large a. However, Rf as t also increases with a. 
Larger satellites migrate more rapidly at larger a than at smaller a. For a larger satellite 
with R = 8 km, the migration rate slows considerably from a = 60 Rp (roughly 20 km yr _1 ) 
to a = 40 Rp (roughly 10 km yr _1 ) to a = 20 Rp (< 1 km yr _1 ). The relative change in the 
migration rate follows theoretical predictions (eqs. [32H33] ). 

The presence of a central binary changes these results modestly (Fig. [11] right columns). 
The most obvious impact of the Pluto-Charon binary is the increased velocity dispersion of 
disk particles. Larger velocity dispersion reduces the effectiveness of torque exchange between 
the satellite and disk particles in the satellite's corotation zone, making it more difficult for 
a satellite to open up a gap in the disk (eq. [30]). Thus, smaller satellites migrate more 
slowly around a binary. There is some evidence from the simulations that larger satellites 
may migrate more rapidly around a binary. In these cases, larger satellites first clear their 
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corotation zone and halt their radial drift. As the simulation proceeds, stirring from the 
binary fee ds some material into the c orotation zone, commencing a kind of "attenuated" 



migration (IBromley fc Kenyonl l2011bl ). We speculate that this behavior may allow larger 



objects to drift more quickly through a circumbinary disk. 

Simulations with a variety of disk masses confirm these general trends. Under the right 
conditions in disks with masses larger than ~ 10 18 g, satellites with R pa 1-20 km can undergo 
fast, type II, or attenuated migration with rates ranging up to 20 km yr _1 . In lower mass 
disks, small satellites cannot clear the gap required to initiate fast migration. Although large 
satellites can form gaps in disks with Md < 10 18 g, the minimum radius required to open 
a gap is larger than the maximum radius for fast migration. Thus, large satellites cannot 
migrate in the fast mode. Instead, these objects migrate a factor of ten more slowly in the 
type II mode. 

To test whether the disk can circularize the orbits of small satellites, we derive predicted 
rates for 4-10 km objects in disks with Md = 1 — 10 x 10 17 g surrounding a single central 
object. For satellites in this size range, damping rates are independent of radius. The 
damping time scales with the disk mass: 

t damp « - « 10 5 (^ry 3 ) yr . (48) 



Even in a low mass disk, damping times are comparable to the formation time for 4-10 km 
satellites. 

Estimating damping rates for satellites orbiting a central binary are complicated by 
precession, orbital resonances, and other dynamical issues. For the small rates indicated by 
simulations with a single central object, accurate estimates require a substantial investment 
of cpu time. When the Pluto-Charon binary has a circular orbit, dynamical friction between 
the satellite and disk particles will still circularize the satellite's orbit. In a non-circular 
Pluto-Charon binary, the damping time probably depends on the time scale for the disk and 
tides to circularize the binary. Because these processes act on similar time scales, deriving 
the circularization time for a satellite orbit is very cpu intensive and beyond the scope of 
this study. 



3.3.2. Migration of Multiple Satellites 

To investigate whether ensembles of proto-satellites can migrate, we examine a second 
suite of simulations (Figure IT2|) . Here, we maintain the same disk surface density distribution 
and focus on satellites with R = 4 km or 10 km orbiting within a wide annulus spanning 
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orbital distances of 35 Rp to 55 Rp. Satellites and disk material orbit a compact Pluto- 
Charon binary with orbital separation a = 5 Rp. For material at 35-55 Rp, stirring by the 
binary has a smaller impact than at 20 Rp. 

The top panel of Figure [12] follows a simulation of five 10 km satellites initially evenly 
spaced in orbital separation. With typical Hill radii of 0.2 Rp, these satellites interact 
weakly among themselves. All have radii larger than Rf as t, migration rates are slow. Modest 
scattering of disk particles yields some inward or outward motion and occasional stronger 
gravitational interactions among the satellites. However, these satellites are fairly stable over 
1000 yr and longer timescales. 

Simulations with ensembles of five 4 km objects produce more complicated outcomes. 
With much smaller Hill radii, these satellites interact very weakly among themselves. At 
the onset of each simulation, however, all begin to migrate in the fast mode. Some migrate 
inward; others migrate outward. When neighboring satellites migrate in the same direction, 
migration is limited by the stirred up wa ke of disk particles left beh ind by the first satellite 



to pass through that portion of the disk ( iBromley fc Kenyonll2011bl ). Satellites do not cross 



these wakes. Thus, groups of small satellites migrating in the same direction drift radially 
inward/outward a small distance before the migration stalls. 

Small satellites migrating in opposite directions lead to more interesting outcomes. Of- 
ten, these satellites scatter one another, sometimes exchanging orbits (Fig.[T2J middle panel). 
Depending on the sense of migration and the state of the disk after scattering, satellites may 
undergo multiple scattering events before settling down in roughly stable configurations 
within a stirred up disk. Sometimes, satellites merge. Because merged satellites are massive, 
they migrate more slowly and end up in a stable configuration. In our suite of simulations, 
mergers are less common (~ 10%) than scattering events. 

Simulations with a mix of 4 km and 10 km satellites also lead to interesting outcomes 
(Fig. [12] lower panel). If a small satellite undergoing fast migration lies between two larger 
satellites, its larger migration rate guarantees that it will catch up to one of its slowly moving 
neighbors. When it does catch up, it may 'bounce off' the wake of its neighbor and begin 
migrating in the opposite direction. In roughly half of all interactions, the larger satellite 
either scatters or merges with the smaller satellite. If scattering places the small satellite on 
an orbit with a pericenter inside a m i n , the central binary can eject the small satellite out of 
the system. 
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3.3.3. Summary of Migration Calculations 

For satellites with R ~ 4-10 km in a disk surrounding the Pluto-Charon system, mi- 
gration is ubiquitous. In a cold, geometrically thin disk, satellites in this size range can 
undergo fast mode or slower, type II mode migration at rates very close to those predicted 
by theory. These rates lead to significant radial motion, ~ 1-10 Rp in 10 3 yr. Scaling these 
rates to the hotter disks expected from the coagulation and hybrid calculations in §3.1-3.2, 
satellites with radii similar to P4 and P5 may migrate significantly as they grow. Although 
larger satellites such as Nix and Hydra migrate more slowly, radial motion of a few Pluto 
radii seems likely after they reach their final sizes. 

Interactions between the disk and the central binary can augment or reduce migration 
rates. Although changes to the rates are modest, smaller (larger) satellites generally migrate 
more slowly (rapidly). Changes are more significant closer to the binary, where the jitter of 
particle orbits in the disk change the number of particles in the corotation zone of a satellite. 
Thus, small (large) satellites which form closer to the binary are less (more) likely to migrate 
than satellites which form farther away from the binary. 

Ensembles of migrating satellites produce a variety of interesting dynamical phenomena. 
In systems with several small satellites or a mix of large and small satellites, differential 
migration enables large-scale scattering events and satellite mergers. When scattering places 
a satellite on a orbit with a pericenter smaller than a m i n , the central binary ejects the satellite 
from the system. Mergers among migrating satellites allow lower mass disks to produce more 
massive satellites. 



4. DISCUSSION 

In §2-3, we have developed the first complete, numerical model for the formation of 
Pluto's small satellites. Our approach examines two plausible origins for the icy building- 
blocks of these satellites, debris from the giant impact that p roduces the P luto-Charon 



binary and material accreted from the protostellar debris disk. ICanupl ( 120111 ) summarizes 
a set of giant impact calculations which yield enough debris to produce the satellites. Our 
analysis indicates that a circumplanetary or circumbinary disk can capture enough mass 
from the protoplanetary disk to enable formation of low mass satellites. If Pluto's satellites 
are more massive than their minimum allowed masses, satellite formation in debris from the 
giant impact is the more likely mechanism. 



In each scenario, solids end up in a disk or ring surrounding the Pluto-Charon binary. 
For the disk/ring masses required to match the total mass in the known satellites, solids 
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rapidly grow into a large ensemble of proto-satellites. Dynamical interactions among the 
proto-satellites begin an extended period of chaotic growth where multiple satellite-satellite 
mergers lead to a stable system consisting of a few satellites. The number of stable satellites 
is inversely correlated with the starting mass: disks/rings with larger (smaller) initial masses 
produce fewer (more) satellites with larger (smaller) masses. 

The timescale for satellite growth depends on the origin of the circumbinary debris 
disk. In a disk with mass and surface density £ oc a -1 , the formation time for 1-10 km 
satellites is 



If satellites form in captured material, most of the mass lies at large distances, a ~ 10 3 Rp. 
Formation time are then comparable to the time scale for tidal forces to circularize and to 
expand the central binary, tf ~ 1 Myr. Satellites growing out of debris from the giant impact 
probably form much closer to the central binary, with tf ~ 10 3 10 4 yr. 

During chaotic growth, proto-satellites radially migrate through the leftover debris. 
Although proto-satellites cannot fall into the Pluto-Charon binary, satellites can traverse 
a substantial fraction of the circumbinary disk throughout chaotic growth. Migration can 
induce mergers among lower mass satellites, enabling more massive satellites to form in lower 
mass disks. Ensembles of satellites with similar masses migrate in step with one another, 
enabling satellites to find stable equilibrium orbits around Pluto-Charon. 

Migration rates depend on the nature of the central object. Stirring by a central binary 
can heat the particle disk and slow migration of satellites with radii of 0(1) km. The 
binary may also increase migration rates of larger satellites. A possible explanation is that 
migration comes from the asymmetry of inward and outward satellite-disk scattering events. 
The binary will preferentially disturb the orbits of the inwardly scattered disk particles and 
prevent their return to the satellite in the same way as in the case of a single central mass. 
We plan additional simulations to explore this possibility. 

These results suggest that the formation of low mass satellites around Pluto-Charon 
is a natural outcome of planet formation. Raw material for satellite production is readily 
available from a giant impact or the protostellar debris disk. Once the raw material collects in 
a disk or ring, satellite formation is ubiquitous. Stable satellites have orbital properties and 
a range of masses reasonably consistent with known satellites surrounding Pluto-Charon. 
Chaotic growth often leads to satellites with period ratios similar to those of the known 
satellites. Migration can also produce stable satellites in approximate resonances. 




(49) 



Together with previous results from 



Youdin et al 



( 120121 ) . our calculations indicate sev- 
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eral conclusions about each satellite. 



The semimajor axis of P5 is almost as close as possible to the minimum stable semima- 
jor axis a m i n . In our framework, P5 has a low mass because (i) the debris ring had less 
mass close to Pluto-Charon (ii) it was scattered out of a more massive region of the 
ring by Nix or Hydra, or (iii) after formation at large a, it migrated into an innermost 
stable orbit. 

The or bit of P4 lies withi n a small stable region in between the orbits of Nix and 
Hydra (lYoudin et al.ll2012l ). More massive satellites in this region are unstable. From 
the coagulation and migration calculations, P4 is a smaller remnant of the accretion 
process. 

The masses of Nix and Hydra are consistent with the coagulation models. The central 
cores of either satellite could have been formed the giant impact and then accreted 
debris leftover from the impact. 



Aside from these basic successes, our framework has one major challenge. 

Once debris is in stable orbits around Pluto-Charon, th e binary excites rad ial and 
epicyclic motions in orbiting particles (e.g., Lee fc Peak 2006 : Youdin et elI.I I2OI2I ) . Both 
motions increase the radial velocity dispersion used in our calculations. For stable orbits 
around a circular Pluto-Charon binary, the radial oscillation is small compared to the velocity 
dispersion. Because collisional damping is so effective at reducing the velocity dispersion, 
oscillations induced by a circular binary are generally unimportant. 

If the central binary has e > 0.01, the epicyclic motion is comparable to the typi- 
cal velocity dispersion during runaway and oligarchic growth. Thus, an initially eccentric 
Pluto-Charon binary could maintain a large velocity dispersion among debris particles, pro- 
mote an extended phase of fragment ation, and frustrate the gr owth of proto-satellites in a 
circumbinary ring (see, for example, iPaardekooper et al.l 12012k in the context of circumbi- 
nary planets). In our calculations, however, collisional damping is very effective among small 
debris particles and collision fragments. Thus, damping and fragmentation can counter the 
impact of a central eccentric binary. 

The long-term eccentricity evolution of the central binary is the main uncertainty in 
evaluating the relative impact of collisional damping and binary stirring. In the s tandard 
tidal picture, e damps on the tidal time scale of 0.1-10 Myr (jDobrovolskis et al.lll997l ). When 
damping is so slow, collisional grinding in a low mass disk surrounding an eccentric binary 
probably reduces the mass in solids considerably. Satellite formation is challenging. In both 
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modes of satellite formation, however, tidal damping of e is augmented by interactions with 
the surrounding particle disk/ring. When a ring lying close to the binary has a mass of 
~ 10 23 g - roughly the maximum inferred from giant impact calculations - the damping 
time scale of ~ 4 x 10 3 yr is much smaller than the tidal time scale. On this short time 
scale, collisional grinding reduces the disk mass to levels similar to the total masses of the 
small satellites. After the binary circularizes, collisional damping reduces particle velocities, 
enabling satellite growth. 

Because rings formed during the giant impact have more mass closer to the central 
binary than disks generated from the protoplan etary disk, rings damp the central binary 



more rapidly (e.g.. lLin &: Papaloizoulll979l . ll986al ). However, satellites growing in a captured 



disk grow too slowly and are too far away from the central binary to care about the binary 
eccentricity or the damping time. Thus, both mechanisms appear viable if the Pluto-Charon 
binary has a large, initial eccentricity. 

Constraining the imp act of the central binary requires a new set of simulations (see also 



Paardekooper et al.ll2012l ). Initial tests are promising and suggest that satellite formation 
is likely even in an eccentric binary system. Aside from these new simulations, other more 
comprehensive numerical simulations can help constrain satellite formation models. 

• In the protostellar accretion and giant impact models, improved SPH simulations of col- 
lision outcomes would clarify the likely range of initial conditions available for satellite 
formation. In giant impact models, better mass resolution would yield the distribution 
of small particle masses within the debris. This result would enable a better evaluation 
of the possibility that small proto-satellites survive the giant impact more or less intact 
and serve as the seeds for Pluto's small satellites. 



In accretion models, more robust SPH models (as in ICanupl (120111 )) for the outcomes 
of (i) two objects colliding within the Hill sphere of a planet and (ii) a single object 
from the protoplanetary disk colliding with an object in the circumplanetary disk would 
provide better constraints on our estimates for the maximum mass of a circumplanetary 
disk and its evolution. These improved estimates would allow us to make a better 
evaluation of the viability of this model for satellite formation. 

In the aftermath of a giant impact, our estimates assume efficient collisional damping 
among small objects in the debris. Improving on this approach requires a detailed 
analytic model or a complete numerical model which shows how the properties of 
the circumbinary ring depend on the initial conditions of the debris and the physical 
interactions between debris particles. 
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Once the impact debris lies in a ring surrounding the binary, many physical processes 
occur on roughly the same time scale. In §3, we consider three separate calculations for 
(i) the expansion of the ring into a disk, (ii) the growth of debris particles into satellites, 
and (iii) the migration of satellites through the leftover disk. Our results demonstrate 
that particles grow as the ring diffuses; proto-satellites migrate through the expanding 
disk of smaller particles as the satellites grow. Combining these three calculations into 
a single calculation requires merging the coagulation, diffusion, and migration codes 
within Orchestra. Our initial tests of this new algorithm are promising. We plan to 
describe a complete suite of simulations with this algorithm in subsequent papers. 

Aside from combining coagulation, diffusion, and migration into a single calculation, 
it is also necessary to include the long-term evolution of the central binary. Detailed 
tidal models show that t he time scale for tides to circularize and expand the binary 
(jDobrovolskis et al .11 1 9 9 71 ) are longer than the diffusion and growth times, but similar 
to the time scale for chaotic growth to merge several dozen proto-satellites into the 
satellites we see today. These time scales are 10 — 10 3 times longer than the time scale 
for a ring of debris to damp the eccentricity and compress the orbit of Pluto-Charon. 
Based on our migration calculations, the interplay of tides, dynamical interactions 
between the debris and Pluto-Charon, and migration likely establishes the current 
orbits of the satellites. Thus, treating these physical processes accurately is essential 
for understanding the current orbital arrangement of the satellites. 



Aside from informing the next generation of numerical simulations for satellite forma- 
tion, the New Horizons mission can test several predictions of our calculations. 



In our simulations, forming four (or more) satellites orbiting Pluto-Charon requires 
relatively small disk masses, < 3 x 10 20 g. Together with numerical simulations of 
the orbital stability of the satellites (lYoudin et al .1120121 ). these results require satellite 
albedos, A 0.4-1. This range for A is similar to the measured A ?s 0.5 for Pluto and 



Charon fsee Marcialis et al.lll992l ; iRoush et al.lll996l : lBrowJ200a iBuie et al.ll2010al lb. 
and references therein) and much larger than t he A for many Kuiper belt objects 
(e.g., iStansberry et al.ll2008l ; iBrucker et al.l 120091 ). Direct measurements of A by New 



Horizons will test these constraints. 

All of the hybrid calculations leave several very small satellites orbiting beyond the 
orbit of Hydra. The radii of typical leftovers - R < 1-3 km - imply optical magnitudes 
of 28 or fainter. Tes t ing t his prediction is a challenge for Hubble Space Telescope 
( IShowalter et al.ll201lL 120121 ). but satellites this size or smaller are easily visible during 
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the flyby of New Horizons. For objects with R m 1-3 km, the most likely range of 
semimajor axes is a ~ 70-90 Rp. Smaller satellites might have a w 100-200 Rp. Both 
sets of small satellites should have small inclinations relative to the Pluto-Charon 
orbital plane, similar to the observed i for the known satellites. If smaller particles in 
a leftover debris disk (see next item) have sufficient mass, these outer satellites should 
also have circular orbits. Smaller masses in leftovers imply larger e for small outer 
satellites. We plan numerical calculations to test the longevity of these satellites in the 
tidal field of the Sun-Neptune system. 

Our calculations predict little or no mass between P5 at a ~ 39 Rp and the innermost 
stable orbit at a m i n ~ 34 Rp (see also lYoudin et al.ll2012h . If this region contains debris 
from the giant impact or from occasional impacts of small Kuiper belt objects with 
the satellites, New Horizons might detect a low density stream of material from a min 
to Charon's orbit at ac ~ 17 Rp. As with circumbinary gaseous disks, the properties 
of this stream might yi eld constraints on the intern al viscosity of small solid objects 
orbiting Pluto-Charon ( lArtymowicz fc Lubowlll996l ). 



All of the diffusion calculations predict ensembles of smaller particles, R < 0.01-10 m, 
in an extended disk well beyond the orbits of the known satellites. For a mass of 
~ 10 15 g in leftover small particles with typical radius R, the predicted optical depth 
of the disk is t\ ~ 10~ 10 R~ 1 at semimajor axes a ~ 70-200 Rp. Larger (smaller) disk 
masses imply proportionately larger (smaller) optical depths. Although this mass is 
significant, the predicted optical depth is much smaller than the r < 10~ 7 estimated 
for clouds of 1-1 00 jum particles produc ed by impacts of Kuiper belt objects on the 
known satellites (IPoppe fc Horanyill201ll ). New Horizons measurements of the optical 
depths of small particles surrounding the known satellites and at larger distances should 
place interesting constraints on the long-term evolution of the satellite system and its 
interactions with the Kuiper belt. 

Measuring the outer radius of the extended disk tests the capture hypothesis. In this 
picture, most of the mass lies at large semimajor axes, a ~ 1000 Rp. Identifying any 
material at a m 100-1000 Rp suggests that capture is a viable mechanism. The amount 
of material at these distances provides a measure of the capture efficiency (see eq. [8]). 

Composition data from New Horizons should also enable tests of the formation scenar- 
ios. In the giant impact picture, satellites form out of the same material as Charon. De 
pend i ng on relative contamin ation from impacting Kuiper belt objects (e.g. lStern et al. 
20061 ; IPoppe &: Hor anvil 1201 ll ). satellites should have roughly the same composition as 
Charon. In the circumplanetary accretion scenario, satellites form out of the collisional 
debris of Kuiper belt objects captured by Pluto, Pluto-Charon, or both. The broad 
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variety of physical characteristics of Kuiper belt objects (see lBarucci et al.l 120081 ) sug- 
gests that satellites grown from the debris of Kuiper belt objects should have very 
different co mpositions from satellites grown from the debris of a giant impact (see also 
Stern! 120091 ). Thus, direct measurement of the compositions of individual satellites and 
gradients in the composition among any smaller satellites and debris should constrain 
both scenarios for satellite formation. 

Observed shapes of satellites might test formation models. In general, large fragments 
from the impact should be more irregular than satellites grown from a circumbinary 
disk/ring of much smaller particles. 



To develop predicted configurations for the Pluto-Charon system during a New Horizons 
encounter, we ran n-body simulations of Pluto-Charon, the known small satellites, and 
several predicted satellites embedded within a disk of two million massless tracer particles. 
After 20 yr (roughly 1000 orbits of the inner binary), the satellites have cleared many tracers 
from the ir orbits fFig. IT5]). On longer time scales, the satellites clear their orbits completely 
(see also lYoudin et al.ll20121 ). The predicted optical depth and extent of the disk depends on 
(i) the number of small satellites beyond Hydra, (ii) the viscosity of disk particles, and (iii) 
the amount of mass in small partic les produced from recent impacts on the satellites (e.g., 



Sterol |2009| ; iPoppe fc Horanyi 



201ll ). If (i) there are few small satellites, (ii) disk particles 



have low viscosity, and (iii) impacts are common, we expect New Horizons to detect a 
modest disk of small particles with rings and gaps similar to those shown in the Fig. If, 
however, (i) there are many small satellites beyond Hydra's orbit, (ii) disk particles have 
large viscosity, and (iii) impacts are relatively rare, the New Horizons encounter may reveal 
little disk material. Once New Horizons has mapped the Pluto-Charon system completely, 
improved numerical simulations will enable a more complete examination of the formation 
and history of the satellite system. 



SUMMARY 



We describe a theoretical framework for the origin of Pluto's small satellites in the 
context of t wo scenarios, (i) a giant impact that produces the Pluto-Char on binary (e.g., 
Canudl201ll ) or (ii) accretion from material in the protoplanetary disk (e.g.. iLithwick fc Wu 
20081). Our main results follow. 



Throughout the early history of the solar system, giant impacts capable of producing 
the Pluto-Charon binary have a frequency of 1 per 100-300 Myr. These impacts 
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leave behind enough debris to enable satellite formation (lCanupll201ll ). After impact, 
Charon captures ~ 50% of the debris. Collisional damping enables debris particles 
ejected by Charon to attain bound orbits at a few times the orbital separation of the 
Pluto-Charon binary. 

Capturing debris from collisions of Kuiper belt objects within the Hill sphere of Pluto 
(prior to impact) or Pluto-Charon (after impact) can produce a circumplanetary disk 
with a maximum mass of roughly 3 x 10 19 g. Unless the growth of satellites in the 
captured material is nearly 100% efficient and the current satellites have albedo Am 1, 
satellites are more likely to form in debris from the impact than in material captured 
from the protoplanetary disk. 

In a ring of small particles, R < 0.1 km, surrounding Pluto-Charon, satellites with 
R 10-80 km form rapidly on a time scale of 10 3 — 10 4 yr. The formation time is 
similar to the migration time scale and the viscous time scale. Thus, satellites grow 
and migrate in a viscously spreading disk surrounding Pluto-Charon. 

After the impact, objects with R m 1-10 km might survive collisional disruption during 
the ejection process. Calculations of disks with small particles and a few seeds produce 
fewer satellites with larger masses than calculations with small particles alone. 

Radial migration encourages mergers of small satellites within an evolving disk sur- 
rounding Pluto-Charon. 

Scattering events among differentially migrating satellites sometimes places objects 
on orbits with pericenters close to the central binary. Unless interactions with disk 
particles can damp the orbit, the Pluto-Charon binary ejects these satellites from the 
system. 

For identical masses of the central planet (s), the surrounding disk, and a small satellite 
embedded in the disk, migration rates around a binary planet differ slightly from those 
around a single planet. For Pluto-Charon, the drift rate of a satellite the size of P5 
maybe be as high as 0.5 km yr -1 . 

For a disk with mass sufficient to produce the satellites, the time scale for dynamical 
interactions to damp e and reduce a for the central binary is shorter than the tidal 
time scale. Thus, an accurate model for the tidal evolution of Pluto-Charon must 
include the evolution of the circumplanetary ring that includes growing and migrating 
satellites. 
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The cur rent number of sma ll Pluto-Charon satellites strongly favors low mass disks 
(see also lYoudin et al.ll2012l ). In our calculations, disks with masses ~ 3 — 10 x 10 19 



have the best chance at producing 3-5 satellites with orbital semimajor axes a w 40- 
60 Rp. Matching observations requires Pluto's satellites to have albedos A m 0.4-1. 
New Horizons data will test this prediction. 

Model satellites produced in formation calculations have orbital i close to those ob- 
served, but the typical e is generally too large. Orbital migration can reduce e to 
acceptable levels. 

The calculations predict a few very small, R < 1-3 km, satellites and an extended disk 
of even smaller particles, R ~ 1-100 cm, beyond the current orbit of Hydra. Detecting 
these satellites and the disk from the ground is very challenging. If they are present, 
New Horizons should detect them easily. 

Material lying between P5 and the innermost stable orbit around Pluto-Charon might 
produce a stream of small particles from the circumbinary disk to Charon's orbit. 



Our results demonstrate that numerical calculations can produce simulated satellite 
systems with properties reasonably close to those observed. Prior to the New Horizons 
flyby of Pluto-Charon, we expect to refine the predictions considerably. Once New Horizons 
probes the masses, orbital architecture, and composition of the Pluto-Charon binary, a rich 
interplay between the data and the numerical simulations will enable a much more robust 
theory for the formation of satellites (planets) in binary planetary (stellar) systems. 



We acknowledge generous allotments of computer time on the NASA 'discover' cluster 
and the SI 'hydra' cluster. Advice and comments from M. Geller and A. Youdin greatly im- 
proved our presentation. Portions of this project were supported by the NASA Astrophysics 
Theory and Origins of Solar Systems programs through grant NNX10AF35G and the NASA 
Outer Planets Program through grant NNX11AM37G. 
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Table 1. Key for Frequently Used Symbols 



Symbol Definition 



R, M, p, A, v esc radius, mass, mass density, albedo, and escape velocity of a solid particle 

a, e, i, fi, T semimajor axis, eccentricity, inclination, angular velocity, and period of an orbit 
v orbital velocity 

Rh Hill radius 

Rp, Mp, Mc, qpc radius, mass of Pluto; mass of Charon; mass ratio of Pluto-Charon 
ape, epc, QpCi Tpc orbital parameters of Pluto-Charon binary 

u), t r angular rotational velocity and rotational period of Pluto or Charon 

t exp , tsync time for Pluto-Charon orbit to expand or to synchronize with rotational period 

S,x m surface density and mass scaling factor, S 

v,t v viscosity, viscous time scale; t v w a 2 jv 

a m i n minimum stable orbital distance from the Pluto-Charon binary 

H z , v z , v r vertical scale height and velocity dispersion; radial velocity dispersion 

Td,T p optical depth of circumplanetary ('d') or protoplanetary ('p') disk; r oc T,/R 

ad, Md, S radial distance, mass, and fractional width of debris ring 

tdf time for dynamical friction to damp orbital eccentricity of disk particles 

R m in, Rmax si zc of smallest, largest particle 

t c , a, f g collision time, collisional cross-section, gravitational focusing factor 

b, Vi mp impact parameter and impact velocity 

p, f probability for giant impact, efficiency factor for accretion 

Q c , Q* D collision energy; energy required to disperse half the mass to infinity 

m pro promotion mass in hybrid code; the n-body code tracks objects with m > m pro 

t P ro,tf time scale for satellites to reach the promotion mass or their final mass 

N s , Np number of satellites or Pluto- mass objects formed in simulations 

Rgap Hill radius of satellite able to open a gap in the disk 

Rfast maximum Hill radius of satellite able to migrate in the fast mode 

a fast , an migration rate in the fast or type II mode 
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Table 2. Physical Processes Involved in Satellite Formation 



Phase 



Physical Process 



Section 



Giant Impact 



Debris Ejection 



Ring Formation 



Satellite Growth 



Migration 



Model overview 

Pluto-Charon mass objects grow in protoplanetary disk 2.1 

Pluto-Charon mass objects collide in disk 2.1 

Circumplanetary disk evolution 2.2 

Pluto-Charon accrete some debris 2.3 

Pluto-Charon eject some debris 2.3 

Collisions circularize orbits of debris particles 2.3 

Tidal forces circularize and expand the binary orbit 2.4 

Debris circularizes binary orbit 2.4 

Viscosity spreads ring into disk 2.5.1 

Collisions destroy largest debris particles 2.5.2 

Collisional damping reduces velocity dispersions 2.5.3 

Satellites begin to grow 2.5.3 

Fast migration of small satellites 2.6 

Slow migration of large satellites 2.6 



Viscous Disk 



Satellite Growth 



Migration 



Numerical Simulations 

Evolution of material captured from protoplanetary disk 3.1.2 

Evolution of debris from the giant impact 3.1.3 

Pure coagulation calculations of small satellites 3.2.1 

Hybrid calculations of large satellites 3.2.2 

Migration of single satellites 3.3.1 

Migration of multiple satellites 3.3.2 

Mergers 3.3.2 
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Table 3. Fragmentation parameters for coagulation calculations 



Type 


Qb 


Pb 


Q 9 


Pg 


Weak 


2 x 10 5 


-0.4 


0.22 


1.30 


Strong 


lO^lO 5 


+0.0 


1.5 


1.25 


Very Strong 


lC^-lO 5 


+0.0 


10~ 4 


2.00 
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Fig. 1.— Probability of giant impacts in protoplanetary disks. Giant impacts require 
collisions between a massive target, m t > 10 25 g, and a projectile with m p > 0.3m t . In each 
panel, the legend lists the slope of the radial surface density profile n and the relative disk 
mass x m . Colored curves illustrate the time evolution of the impact probability per Myr for 
disk annuli at a = 30-36 AU (violet), a = 36-44 AU (blue), a = 44-54 AU (turquoise), a 
= 54-66 AU (green), a = 66-80 AU (orange), a = 80-100 AU (magenta), a = 100-120 AU 
(maroon), and a = 120-150 AU (black). Giant impacts happen earlier and more frequently 
in more massive disks. 
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Fig. 2. — As in Fig. [T] for disks with n = 1.5 and x m = 1.00. Left panels: results for 
calculations with Rq = 1 km and different fragmentation parameters. Right panels: results 
for calculations with weak planetesimals and Rq = 1, 10, and 100 km. Because larger objects 
form more frequently when planetesimals are strong, giant impacts are more common in 
calculations with strong planetesimals. Planets form faster when planetesimals are initially 
smaller; thus, giant impacts occur earlier when Rq is smaller. 
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Fig. 3.— Surface density profiles in a circumplanetary disk surrounding a Pluto- mass 
protoplanet. The disk has a gaussian initial surface density with a mass of 10 11 g centered at 
20 Rp. Collisions from material in the protoplanetary disk augment and deplete the surface 
density according to the prescription in eq. f )39|) . Mass addition gradually erases the initial 
conditions, yielding a roughly constant surface density at all radii at 10 4 yr (violet curve), 
10 5 yr (blue curve), 10 6 yr (green curve), 10 7 yr (orange curve), and 10 8 yr (magenta curve). 
As the disk mass reaches a steady-state at 10 9 yr (maroon curve), the surface density peaks 
close to the central protoplanet. 
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Fig. 4. — Surface density evolution in a circumplanetary disk with = 3 x 10 20 g sur- 
rounding the Pluto-Charon binary at 10 4 yr (upper left panel), 10 5 yr (upper right panel), 
10 6 yr (lower left panel), and 10 7 yr (lower right panel). The calculations assume a constant 
ape = 5 Rp (vertical dashed line in each panel). The vertical grey bar indicates separations 
occupied by the satellites P5, Nix, P4, and Hydra. Solid colored lines indicate the surface 
density distribution for 0.1 km particles with radial velocity dispersion v r = 5 cm s _1 (blue), 
50 cm s _1 (orange), and 500 cm s _1 (maroon). The binary truncates the disk at a distance 
of 2 ape (10 Rp) from Pluto. Disks with large v r expand more rapidly than disks with small 
v r . The outer radii of the disks reach the orbits of the small satellites in 10 5 -10 7 yr. 
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Fig. 5. — As in Fig. H]for an expanding Pluto-Charon binary. In each panel, the dashed 
vertical line indicates the semimajor axis of Charon's orbit. Compared to calculations with- 
out an expanding binary, the disks in these calculations are more compressed but reach the 
orbits of the small satellites at similar times. 



-65 - 



2 - 



1 ih 

s 

OS 

OX) 

5 Oh 



-1 = 



3 x 10 21 g model 




3 x 10 19 g model 



log Time (yr) 



Fig. 6. — Radius evolution for satellites in a disk surrounding the Pluto-Charon binary. 
The dashed lines indicate the lower limit for the radius of P5 {A — 1) and the upper limit 
for the radius of Hydra (A = 0.04) assuming a mass density p = 1 g cm -3 . The thin (thick) 
lines plot results for a model with an initial disk mass of Md = 3 x 10 19 g {Md = 3 x 10 21 g) 
in solid objects with initial radii of 0.1 km and smaller. Each line indicates the radius 
evolution for the largest object in annuli at 17.5 Rp (violet), 20 Rp (blue), 23 Rp (cyan), 
26.5 Rp (turquoise), 30.5 Rp (green), 35 Rp (orange), 40.5 Rp (maroon), and 46.5 Rp 
(brown). Satellites grow faster in more massive disks. The range of model radii span the 
range inferred for satellites of Pluto-Charon. 
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Fig. 7. — Semimajor axis evolution for satellites in a disk surrounding Pluto-Charon. The 
upper left corner of each panel lists the initial disk mass. Along the right edge of each panel, 
numbers indicate the radii (in km) of each satellite. More massive disks produce fewer, larger 
satellites than less massive disks. 
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Fig. 8. — As in Fig. [TJfor an initial disk mass ~ 3 x 10 19 g composed of 0.1 km and smaller 
planetesimals. Upper panel: evolution without larger planetesimals; lower panel: evolution 
with two 10 km planetesimals. Calculations with a several initial large planetesimals yield 
fewer, larger satellites than calculations without large planetesimals. 
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Fig. 9. — Radius evolution for n-bodies. The calculation begins with rrid = 3 x 10 21 g and 
10 18 g. Once a few large objects form, they gradually accrete all of the smaller 

n-bodies. 
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Fig. 10. — Distribution of orbital eccentricity e and inclination i for satellites produced in 
disks with initial masses of 3 x 10 19 g (violet points), 3 x 10 20 g (blue points), and 3 x 10 21 g 
(orange points). Although there is no clear correlation of average orbital inclination with 
disk mass, disks with larger initial masses produce satellites with larger average e. The 
horizontal and vertical bars indicate the average of the measured eccentricity and inclination 
for Nix (e^r and zjy) and Hydra (e# and z#) from the 2007 Jacobson PLU017 JPL satellite 
ephemeris. For the model satellites, 10% have e < 6h,n an d 83% have i < ih,n- 



-70 - 



Q 



CO 

a: 




no binary, 

= 20 R p 
40 R p 
60 R p 




a re = 5 R P' r P hy = 8 km 



a = 20 R p 
40 R p 
60 R p 




100 200 300 400 
Time (yr) 



100 200 300 400 500 
Time (yr) 



Fig. 11. — Migration of single satellites through a particle disk. Each trace shows the result 
of a simulation with either a single central point mass (left panels, with mass equal to the 
sum of Mp and Mc) or a central Pluto-Charon binary with orbital separation of 5Rp. The 
legends give the physical radii of each satellite with p — 1 g/cm 3 ; the legends also show the 
semimajor axis of the satellites' initial orbit. The disk model has a total mass of 3 x 10 20 g 
in a disk with surface density E ~ 1/a extending in radius from 20Rp to 70Rp. Note that 
migration is more commonly in the inward direction, however, it can be outward by chance. 
Solid curves show inward drift, while dashed lines refer to outward motion. 
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Fig. 12. — Migration of multiple satellites around a tight Pluto-Charon binary. The binary 
orbit and disk model are the same as in Figure [TTJ except that each set of five satellites is 
embedded in a 20Rp annulus of disk particles. Each panel gives an example of the simulation 
output; the upper panel shows equal- mass 10 km satellites, while the middle panel tracks 
smaller 4 km bodies. In the former case, some radial drift occurs but does not lead to 
significant scattering or merging. The smaller satellites are more interactive; the middle 
panel shows an orbit crossing event. The lower panel illustrates differential migration in 
which a smaller body merges with a larger one. 



Fig. 13. — Predicted system configuration for the Pluto-Charon system. The Pluto-Charon 
binary (represented by the two largest white disks), the four small satellites - P5, Nix, 
P4, and Hydra (represented by the four small white disks), and three smaller satellites 
(represented by the green disks) lie within an extended ensemble of solid particles shown as 
small blue dots. This configuration is the result after twenty years of a computer simulation 
with two million massless tracer particles surrounding the known and predicted moons. On 
this short time scale, the small satellites clear out most of the tracers along their orbits. On 
much longer time scales, satellites will clear tracers from a larger fraction of their orbits. 



